Benchmarking clinical risk prediction algorithms with ensemble machine learning: An illustration of the superlearner algorithm for the non-invasive diagnosis of liver fibrosis in non-alcoholic fatty liver disease
Vivek Charu, Jane W. Liang, Ajitha Mannalithara, Allison Kwong, Lu Tian, and W. Ray Kim
2024-03-11
# Load libraries
library(tidyverse)
library(gridExtra)
library(gtable)
library(reshape2)
library(forcats)
library(cowplot)
library(SuperLearner)
library(cvAUC)
library(pROC)
library(survey)
library(DescTools)
library(WeightedROC)
# Set background to be white for all ggplots
theme_set(theme_classic())Data processing
We read in the training data (NASH-CRN) and two validation sets (FLINT and NHANES).
# Data directory
dir = "/Users/jliang/Library/CloudStorage/Box-Box/Jane-Vivek/Superlearner/Data/"
# NASH-CRN
nash = read.csv(paste0(dir, "NASH_CRN/NASH.CRN.table3.csv"))
# FLINT
flint = read.csv(paste0(dir, "FLINT/FLINT_Tbls_1and2.csv"))
# NHANES
nhanes = read.csv(paste0(dir, "NHANES_V2/Nhanes_2017_20_ExamData_SAFE_V2.csv"))The names of the NASH-CRN variables that will be used to fit the superlearner models are defined as follows. We also name the outcome variable and distinguish between continuous variables (which can be log transformed) and binary variables (which will not be log transformed).
# Predictors to use for model specification
pred_vars = c("AST", "ALT", "ALKA", "BILIT", "GGT", "GLUC", "ALB", "HEMA",
"WBC", "PLAT", "CHOL", "TRI", "HDL", "LDL", "HBA1C", "AGE",
"MALE", "HISPANIC", "WHITE", "BMI", "DIAB2", "HTN", "GLOB")
# Continuous predictors to log transform
pred_vars_log = c("AST", "ALT", "ALKA", "BILIT", "GGT", "GLUC", "ALB", "HEMA",
"WBC", "PLAT", "CHOL", "TRI", "HDL", "LDL", "HBA1C", "AGE",
"BMI", "GLOB")
# Binary/categorical predictors that will not be transformed
pred_vars_binary = c("MALE", "HISPANIC", "WHITE", "DIAB2", "HTN")
# Outcome variable
response = c("sigfibro2", "sigfibro3")We format the FLINT variables to match the names and units in the NASH-CRN data.
flint$GLUC = flint$GLUCB
flint$PLAT = flint$PLATELET/1000
flint$TRI = flint$TRIG
flint$HISPANIC = flint$HISP
flint$WHITE = ifelse(flint$RACE == "White", 1, 0)
flint$DIAB2 = flint$DIABRX
flint$GLOB = flint$PROT - flint$ALB We format the NHANES variables to match the names and units in the NASH-CRN data.
nhanes$AST = nhanes$LBXSASSI
nhanes$ALT = nhanes$LBXSATSI
nhanes$ALKA = nhanes$LBXSAPSI
nhanes$BILIT = nhanes$LBDSTBSI/17.1
nhanes$GGT = nhanes$LBXSGTSI
nhanes$GLUC = nhanes$LBXGLU
nhanes$ALB = nhanes$LBDSALSI/10
nhanes$HEMA = nhanes$LBXHCT
nhanes$WBC = nhanes$LBXWBCSI
nhanes$PLAT = nhanes$LBXPLTSI
nhanes$CHOL = nhanes$LBXTC
nhanes$TRI = nhanes$LBXTR
nhanes$HDL = nhanes$LBDHDD
nhanes$LDL = nhanes$LBDLDL
nhanes$HBA1C = nhanes$LBXGH
nhanes$AGE = nhanes$RIDAGEYR
nhanes$MALE = ifelse(nhanes$RIAGENDR == 1, 1, 0)
nhanes$HISPANIC = ifelse(nhanes$RIDRETH3 %in% c(1, 2), 1, 0)
nhanes$WHITE = ifelse(nhanes$RIDRETH3 == 3, 1, 0)
nhanes$BMI = nhanes$BMXBMI
nhanes$DIAB2 = nhanes$Diabetes
nhanes$DIAB2[is.na(nhanes$DIAB2)] = 0
nhanes$HTN = ifelse(nhanes$BPQ020 == 1 | nhanes$BPQ030 == 1 |
nhanes$BPD035 <= 80 | nhanes$BPQ040A == 1 |
nhanes$BPQ050A == 1 | nhanes$BPXOSY1 > 140 |
nhanes$BPXOSY2 > 140 | nhanes$BPXOSY3 > 140 |
nhanes$BPXODI1 > 90 | nhanes$BPXODI2 > 90 |
nhanes$BPXODI3 > 90, 1, 0)
nhanes$HTN[is.na(nhanes$HTN)] = 0
nhanes$GLOB = nhanes$LBXSGB
nhanes$LSM = nhanes$LUXSMED
nhanes$CAP = nhanes$LUXCAPMFor each dataset, we can define a BMI variable that sets values above
40 to 40. The BMI40 variable will be used to calculate the
SAFE score.
# NASH-CRN
nash$BMI40 = nash$BMI
nash$BMI40[nash$BMI40 > 40] = 40
# FLINT
flint$BMI40 = flint$BMI
flint$BMI40[flint$BMI40 > 40] = 40
# NHANES
nhanes$BMI40 = nhanes$BMI
nhanes$BMI40[nhanes$BMI40 > 40] = 40We specify the response variable (significant fibrosis) for each dataset at F2 or higher fibrosis, or F3 or higher fibrosis. For NHANES, this can be defined as those with liver stiffness measure (LSM) above 8 kPa and 12 kPa, respectively.
# Significant fibrosis defined as F2 or higher
nash$sigfibro2 = ifelse(nash$NFIBRO >= 2 , 1, 0)
flint$sigfibro2 = ifelse(flint$BLFIBRO >= 2, 1, 0)
nhanes$sigfibro2 = ifelse(nhanes$LSM > 8, 1, 0)
# Significant fibrosis defined as F3 or higher
nash$sigfibro3 = ifelse(nash$NFIBRO >= 3 , 1, 0)
flint$sigfibro3 = ifelse(flint$BLFIBRO >= 3, 1, 0)
nhanes$sigfibro3 = ifelse(nhanes$LSM > 12, 1, 0)We now drop observations with missing values. This results in a reduction of 693 to 648 observations for NASH-CRN and 283 to 270 for FLINT.
# Subset relevant predictors and drop observations with missing values
# NASH-CRN
nash = na.omit(nash[,c(pred_vars, response, "BMI40")])
# FLINT
flint = na.omit(flint[,c(pred_vars, response, "BMI40")])The NHANES data was collected using a clustered survey sampling
approach. This means that point estimates need to take into account the
sampling weights, and standard errors need to take into account the
clusters and strata. After dropping ineligible NHANES participants using
the Elig_NAFLD variable, 2636 out 14300 observations
remain. Dropping observations with missing values further reduces the
sample size to 1279. Dropping observations with zero sampling weights
results in 1244 observations in the final NHANES-NAFLD dataset.
# NHANES-NAFLD
# Drop ineligible values
nhanes = nhanes %>% subset(Elig_NAFLD == 1)
# Drop observations with missing values
nhanes = na.omit(nhanes[, unique(c(pred_vars, response,
"BMI40", "LSM", "CAP",
"WTSAFPRP", "SDMVPSU", "SDMVSTRA"))])
# Drop people with zero weights
nhanes = nhanes[nhanes$WTSAFPRP!=0,]Some of the superlearners will use log-transformed predictors, so we create log-transformed versions of all continuous variables. One value in the NASH-CRN is infinite, so we replaced it with the mean.
# Log variables
for (var in pred_vars_log) {
nash[,paste0("LN_", var)] = log(nash[,var])
flint[,paste0("LN_", var)] = log(flint[,var])
nhanes[,paste0("LN_", var)] = log(nhanes[,var])
}
# Fix an ALKA value which is zero
nash$LN_ALKA[which(nash$LN_ALKA<(-1000))] = log(mean(nash$ALKA))We define the covariate matrices and response vectors for each dataset. For the NHANES-NAFLD data, we also define the sampling weights, cluster IDs, and strata to use for the complex survey design.
# NASH-CRN
nashX = nash[,pred_vars]
nashX_log = nash[,c(pred_vars_binary, paste0("LN_", pred_vars_log))]
nashX_all = cbind(nashX, nashX_log[,paste0("LN_", pred_vars_log)])
nashY2 = nash$sigfibro2
nashY3 = nash$sigfibro3
# FLINT
flintX = flint[,pred_vars]
flintX_log = flint[,c(pred_vars_binary, paste0("LN_", pred_vars_log))]
flintX_all = cbind(flintX, flintX_log[,paste0("LN_", pred_vars_log)])
flintY2 = flint$sigfibro2
flintY3 = flint$sigfibro3
# NHANES-NAFLD
nhanesX = nhanes[,pred_vars]
nhanesX_log = nhanes[,c(pred_vars_binary, paste0("LN_", pred_vars_log))]
nhanesX_all = cbind(nhanesX, nhanesX_log[,paste0("LN_", pred_vars_log)])
nhanesY2 = nhanes$sigfibro2
nhanesY3 = nhanes$sigfibro3
# Complex survey information and design
nhanes_weights = nhanes$WTSAFPRP
nhanes_clusters = nhanes$SDMVPSU
nhanes_strata = nhanes$SDMVSTRA
nhanes_design = svydesign(id = ~nhanes_clusters,
strata = ~nhanes_strata,
weights = ~nhanes_weights, nest=TRUE)Descriptive statistics
We create a “Table 1” that reports the median and IQR for continuous variables and the mean and standard deviation for binary variables in each of the three datasets. Note that the descriptive statistics for NHANES-NAFLD are weighted by the sampling weights, and the standard deviations also incorporate information from the survey clusters and strata.
# Data dictionary for nice variable names
dat_dict = data.frame(rbind(
c("ALB", "Albumin (g/dL)"),
c("ALT", "Alanine aminotransferase (U/L)"),
c("ALKA", "Alkaline phosphatase (U/L)"),
c("AST", "Aspartate aminotransferase (U/L)"),
c("BMI", "Body mass index (kg/m2)"),
c("GGT", "Gamma-glutamyl transferase (U/L)"),
c("GLOB", "Globulin (g/dL)"),
c("GLUC", "Glucose (mg/dL)"),
c("HDL", "HDL cholesterol (mg/dL)"),
c("HEMA", "Hematocrit (%)"),
c("HBA1C", "Hemoglobin A1C (%)"),
c("HTN", "Hypertension (yes/no)"),
c("LDL", "LDL cholesterol (mg/dL)"),
c("PLAT", "Platelets (1000/mm3)"),
c("BILIT", "Total bilirubin (mg/dL)"),
c("CHOL", "Total cholesterol (mg/dL)"),
c("TRI", "Triglyercides (mg/dL)"),
c("DIAB2", "Type 2 diabetes (yes/no)"),
c("WBC", "White blood cell (1000/mm3)"),
c("AGE", "Age (years)"),
c("HISPANIC", "Hispanic ethnicity (yes/no)"),
c("MALE", "Male gender (yes/no)"),
c("WHITE", "White race (yes/no)"),
c("sigfibro2", "Significant fibrosis of F2 or higher (yes/no)"),
c("sigfibro3", "Significant fibrosis of F3 or higher (yes/no)")
))
names(dat_dict) = c("short", "long")
# Center (median or mean) of each variable
center_tab = t(sapply(c(response, sort(pred_vars)), function(var) {
if (var %in% intersect(dat_dict$short, pred_vars_log)) { # Continuous
c(median(nashX[,var]), median(flintX[,var]),
Median(nhanesX[,var], weights = nhanes_weights))
} else { # Binary variables
if (var %in% response) { # Special case for response variables
c(mean(nash[,var]), mean(flint[,var]),
weighted.mean(nhanes[,var], nhanes_weights))
} else {
c(mean(nashX[,var]), mean(flintX[,var]),
weighted.mean(nhanesX[,var], nhanes_weights))
}
}
}))
colnames(center_tab) = c("NASH-CRN", "FLINT", "NHANES-NAFLD")
# Spread (IQR or standard deviation) of each variable
spread_tab = t(sapply(c(response, sort(pred_vars)), function(var) {
if (var %in% intersect(dat_dict$short, pred_vars_log)) { # Continuous variables
c(IQR(nashX[,var]), IQR(flintX[,var]),
IQRw(nhanesX[,var], weights = nhanes_weights))
} else { # Binary variables
if (var %in% response) { # Special case for response variablee
c(sd(nash[,var]), sd(flint[,var]),
sqrt(svyvar(nhanes[,var], nhanes_design)))
} else {
c(sd(nashX[,var]), sd(flintX[,var]),
sqrt(svyvar(nhanesX[,var], nhanes_design)))
}
}
}))
colnames(spread_tab) = c("NASH-CRN", "FLINT", "NHANES-NAFLD")
# Create Table 1
tab1 = matrix(paste0(round(center_tab, 2), " (", round(spread_tab, 2), ")"),
nrow = nrow(center_tab), ncol = ncol(center_tab))
rownames(tab1) = rownames(center_tab)
colnames(tab1) = colnames(center_tab)
# Use more descriptive variable names
tab1 = merge(dat_dict, tab1, by.x = "short", by.y = 0, sort = FALSE)
rownames(tab1) = tab1$long
tab1$short = NULL; tab1$long = NULL
write.csv(tab1, file = "figs_and_tabs/tab1.csv")
knitr::kable(tab1)| NASH-CRN | FLINT | NHANES-NAFLD | |
|---|---|---|---|
| Albumin (g/dL) | 4.2 (0.5) | 4.3 (0.6) | 4 (0.4) |
| Alanine aminotransferase (U/L) | 64.5 (54) | 68 (61.5) | 21 (15) |
| Alkaline phosphatase (U/L) | 80 (36) | 76.5 (34.75) | 74 (26) |
| Aspartate aminotransferase (U/L) | 46 (34) | 51 (39.5) | 19 (8) |
| Body mass index (kg/m2) | 33.66 (8.6) | 33.58 (7.29) | 32.1 (8.7) |
| Gamma-glutamyl transferase (U/L) | 49 (54) | 48.5 (58) | 24 (16) |
| Globulin (g/dL) | 3 (0.7) | 3 (0.6) | 3 (0.5) |
| Glucose (mg/dL) | 95 (24) | 103.5 (31) | 109 (21) |
| HDL cholesterol (mg/dL) | 42 (14) | 42 (14) | 45 (16) |
| Hematocrit (%) | 42.4 (4.8) | 41.05 (4.88) | 43 (5.3) |
| Hemoglobin A1C (%) | 5.7 (0.9) | 6.2 (1.2) | 5.7 (0.7) |
| Hypertension (yes/no) | 0.44 (0.5) | 0.61 (0.49) | 0.54 (0.5) |
| LDL cholesterol (mg/dL) | 119 (48) | 108.5 (52.75) | 112 (48) |
| Platelets (1000/mm3) | 245 (80.5) | 232.5 (74.5) | 237 (77) |
| Total bilirubin (mg/dL) | 0.7 (0.4) | 0.6 (0.4) | 0.4 (0.3) |
| Total cholesterol (mg/dL) | 193 (51.25) | 185 (63.25) | 183 (52) |
| Triglyercides (mg/dL) | 149 (99) | 154 (92.5) | 116 (79) |
| Type 2 diabetes (yes/no) | 0.22 (0.42) | 0.49 (0.5) | 0.25 (0.43) |
| White blood cell (1000/mm3) | 6.7 (2.4) | 6.85 (2.7) | 7 (2.4) |
| Age (years) | 49 (17) | 53 (16.75) | 54 (26) |
| Hispanic ethnicity (yes/no) | 0.14 (0.35) | 0.16 (0.36) | 0.19 (0.4) |
| Male gender (yes/no) | 0.38 (0.49) | 0.34 (0.48) | 0.55 (0.5) |
| White race (yes/no) | 0.8 (0.4) | 0.82 (0.38) | 0.62 (0.49) |
| Significant fibrosis of F2 or higher (yes/no) | 0.45 (0.5) | 0.59 (0.49) | 0.15 (0.36) |
| Significant fibrosis of F3 or higher (yes/no) | 0.25 (0.44) | 0.33 (0.47) | 0.06 (0.24) |
Train superlearners
For each outcome, we fit superlearners with 12 different base models to all untransformed predictors; all log-transformed predictors; and both untransformed and log-transformed predictors in the NASH-CRN dataset. Default tuning parameters were used for all 12 base models.
# Library of 12 base models
SL_library_12 = c("SL.bayesglm", # Bayesian generalized linear model
"SL.earth", # Multivariate adaptive regression splines
"SL.gam", # Generalized additive model
"SL.gbm", # Generalized boosted model
"SL.glm", # Generalized linear model
"SL.glmnet", # Regularized generalized linear model
"SL.ipredbagg", # Bagging trees
"SL.nnet", # Neural network
"SL.polymars", # Multivariate adaptive polynomial spline regression
"SL.randomForest", # Random forest
"SL.rpart", # Recursive partitioning tree
"SL.svm") # Support vector machine
# F2 outcome
# Untransformed predictors
set.seed(1)
SL_12_F2 = SuperLearner(Y = nashY2, X = nashX, family = binomial(),
SL.library = SL_library_12, method = "method.AUC")
# Log-transformed predictors
set.seed(1)
SL_12_log_F2 = SuperLearner(Y = nashY2, X = nashX_log, family = binomial(),
SL.library = SL_library_12, method = "method.AUC")
# Untransformed and log-transformed predictors
set.seed(1)
SL_12_all_F2 = SuperLearner(Y = nashY2, X = nashX_all, family = binomial(),
SL.library = SL_library_12, method = "method.AUC")
save(SL_12_F2, SL_12_log_F2, SL_12_all_F2,
file = "model_output/SL_12_F2.rData")
# F3 outcome
# Untransformed predictors
set.seed(1)
SL_12_F3 = SuperLearner(Y = nashY3, X = nashX, family = binomial(),
SL.library = SL_library_12, method = "method.AUC")
# Log-transformed predictors
set.seed(1)
SL_12_log_F3 = SuperLearner(Y = nashY3, X = nashX_log, family = binomial(),
SL.library = SL_library_12, method = "method.AUC")
# Untransformed and log-transformed predictors
set.seed(1)
SL_12_all_F3 = SuperLearner(Y = nashY3, X = nashX_all, family = binomial(),
SL.library = SL_library_12, method = "method.AUC")
save(SL_12_F3, SL_12_log_F3, SL_12_all_F3,
file = "model_output/SL_12_F3.rData")We then fit superlearners with 90 different combinations of base models and tuning parameters to the same outcomes and 3 combinations of predictors as above.
# Multivariate adaptive regression spline models (3)
earth_learners =
create.Learner("SL.earth",
tune = list(degree = c(1, 2, 3)))
# Generalized boosted models (27)
gbm_learners =
create.Learner("SL.gbm",
tune = list(interaction.depth = c(1, 2, 3),
shrinkage = c(0.1, 0.01, 0.001),
n.minobsinnode = c(1, 10, 20)))
# Regularized generalized linear models (3)
glmnet_learners =
create.Learner("SL.glmnet",
tune = list(alpha = c(0, 0.5, 1)))
# Bagging trees (9)
ipredbagg_learners =
create.Learner("SL.ipredbagg",
tune = list(minsplit = c(10, 20, 30),
cp = c(0.005, 0.01, 0.05)))
# Neural networks (9)
nnet_learners =
create.Learner("SL.nnet",
tune = list(size = c(2, 5, 10),
decay = c(0, 0.25, 0.5)))
# Multivariate adaptive polynomial spline regression models (9)
polymars_learners =
create.Learner("SL.polymars",
tune = list(gcv = c(2, 4, 6),
knot.space = c(2, 3, 4)))
# Random forests (9)
randomForest_learners =
create.Learner("SL.randomForest",
tune = list(mtry = floor(sqrt(ncol(nashX)) * c(0.5, 1, 2)),
nodesize = c(10, 20, 50)))
# Recursive partitioning trees (9)
rpart_learners =
create.Learner("SL.rpart",
tune = list(minsplit = c(10, 20, 30),
cp = c(0.005, 0.01, 0.05)))
# Support vector machines (9)
svm_learners =
create.Learner("SL.svm",
tune = list(gamma = (1/ncol(nashX)) * c(0.5, 1, 2),
cost = c(0.5, 1, 1.5)))
# Library of 90 base models with varying tuning parameters
SL_library_90 = c(
"SL.bayesglm", # Bayesian generalized linear model (1)
earth_learners$names, # Multivariate adaptive regression spline models (3)
"SL.gam", # Generalized additive model (1)
gbm_learners$names, # Generalized boosted models (27)
"SL.glm", # Generalized linear model (1)
glmnet_learners$names, # Regularized generalized linear models (3)
ipredbagg_learners$names, # Bagging trees (9)
nnet_learners$names, # Neural networks (9)
polymars_learners$names, # Multivariate adaptive polynomial spline regression models (9)
randomForest_learners$names, # Random forest models (9)
rpart_learners$names, # Recursive partitioning trees (9)
svm_learners$names) # Support vector machines (9)
# F2 outcome
# Untransformed predictors
set.seed(1)
SL_90_F2 = SuperLearner(Y = nashY2, X = nashX, family = binomial(),
SL.library = SL_library_90, method = "method.AUC")
# Log-transformed predictors
set.seed(1)
SL_90_log_F2 = SuperLearner(Y = nashY2, X = nashX_log, family = binomial(),
SL.library = SL_library_90, method = "method.AUC")
# Untransformed and log-transformed predictors
set.seed(1)
SL_90_all_F2 = SuperLearner(Y = nashY2, X = nashX_all, family = binomial(),
SL.library = SL_library_90, method = "method.AUC")
save(SL_90_F2, SL_90_log_F2, SL_90_all_F2,
file = "model_output/SL_90_F2.rData")
# F3 outcome
# Untransformed predictors
set.seed(1)
SL_90_F3 = SuperLearner(Y = nashY3, X = nashX, family = binomial(),
SL.library = SL_library_90, method = "method.AUC")
# Log-transformed predictors
set.seed(1)
SL_90_log_F3 = SuperLearner(Y = nashY3, X = nashX_log, family = binomial(),
SL.library = SL_library_90, method = "method.AUC")
# Untransformed and log-transformed predictors
set.seed(1)
SL_90_all_F3 = SuperLearner(Y = nashY3, X = nashX_all, family = binomial(),
SL.library = SL_library_90, method = "method.AUC")
save(SL_90_F3, SL_90_log_F3, SL_90_all_F3,
file = "model_output/SL_90_F3.rData")We summarize the base model coefficients/weights for each superlearner. For brevity, we summed up all weights corresponding to the same model (but with different parameters) in the superlearners with 90 base models.
# Dictionary for nice names of superlearner base algorithms
sl_dict = data.frame(rbind(
c("bayesglm", "Bayesian generalized linear model"),
c("earth", "Multivariate adaptive regression splines"),
c("gam", "Generalized additive model"),
c("gbm", "Generalized boosted model"),
c("glm", "Generalized linear model"),
c("glmnet", "Regularized generalized linear model"),
c("ipredbagg", "Bagging trees"),
c("nnet", "Neural network"),
c("polymars", "Multivariate adaptive polynomial spline regression"),
c("randomForest", "Random forest"),
c("rpart", "Recursive partitioning tree"),
c("svm", "Support vector machine")
))
names(sl_dict) = c("short", "long")
# F2 outcome
# Coefficients for superlearner with 12 base algorithms
coef_SL_12_F2 = cbind("SL-12" = SL_12_F2$coef,
"SL-12 (log)" = SL_12_log_F2$coef,
"SL-12 (all)" = SL_12_all_F2$coef)
rownames(coef_SL_12_F2) = sub("_All", "",
sub("SL.", "", rownames(coef_SL_12_F2)))
# Coefficients for superlearner with 90 base algorithms
coef_SL_90_F2 = cbind("SL-90" = SL_90_F2$coef,
"SL-90 (log)" = SL_90_log_F2$coef,
"SL-90 (all)" = SL_90_all_F2$coef)
rownames(coef_SL_90_F2) = sub("_.*", "",
sub("SL.", "", rownames(coef_SL_90_F2)))
# Combine (sum) weights from same base model type
coef_SL_90_F2 = aggregate(coef_SL_90_F2, list(rownames(coef_SL_90_F2)), sum)
# Put all superlearner coefficients together
coef_SL_F2_tab = merge(coef_SL_12_F2, coef_SL_90_F2,
by.x = 0, by.y = "Group.1")
rownames(coef_SL_F2_tab) = coef_SL_F2_tab$Row.names
coef_SL_F2_tab = coef_SL_F2_tab[,-1]
knitr::kable(round(coef_SL_F2_tab, 3))| SL-12 | SL-12 (log) | SL-12 (all) | SL-90 | SL-90 (log) | SL-90 (all) | |
|---|---|---|---|---|---|---|
| bayesglm | 0.089 | 0.059 | 0.105 | 0.021 | 0.000 | 0.002 |
| earth | 0.246 | 0.045 | 0.007 | 0.064 | 0.000 | 0.006 |
| gam | 0.000 | 0.001 | 0.000 | 0.007 | 0.000 | 0.001 |
| gbm | 0.141 | 0.136 | 0.134 | 0.324 | 0.210 | 0.282 |
| glm | 0.074 | 0.039 | 0.003 | 0.013 | 0.002 | 0.001 |
| glmnet | 0.047 | 0.069 | 0.136 | 0.026 | 0.022 | 0.006 |
| ipredbagg | 0.083 | 0.113 | 0.119 | 0.083 | 0.109 | 0.050 |
| nnet | 0.087 | 0.026 | 0.137 | 0.120 | 0.103 | 0.187 |
| polymars | 0.136 | 0.263 | 0.134 | 0.203 | 0.269 | 0.347 |
| randomForest | 0.095 | 0.146 | 0.129 | 0.034 | 0.160 | 0.078 |
| rpart | 0.003 | 0.002 | 0.001 | 0.007 | 0.011 | 0.017 |
| svm | 0.000 | 0.101 | 0.095 | 0.098 | 0.112 | 0.022 |
# F3 outcome
# Coefficients for superlearner with 12 base algorithms
coef_SL_12_F3 = cbind("SL-12" = SL_12_F3$coef,
"SL-12 (log)" = SL_12_log_F3$coef,
"SL-12 (all)" = SL_12_all_F3$coef)
rownames(coef_SL_12_F3) = sub("_All", "",
sub("SL.", "", rownames(coef_SL_12_F3)))
# Coefficients for superlearner with 90 base algorithms
coef_SL_90_F3 = cbind("SL-90" = SL_90_F3$coef,
"SL-90 (log)" = SL_90_log_F3$coef,
"SL-90 (all)" = SL_90_all_F3$coef)
rownames(coef_SL_90_F3) = sub("_.*", "",
sub("SL.", "", rownames(coef_SL_90_F3)))
# Combine (sum) weights from same base model type
coef_SL_90_F3 = aggregate(coef_SL_90_F3, list(rownames(coef_SL_90_F3)), sum)
# Put all superlearner coefficients together
coef_SL_F3_tab = merge(coef_SL_12_F3, coef_SL_90_F3,
by.x = 0, by.y = "Group.1")
rownames(coef_SL_F3_tab) = coef_SL_F3_tab$Row.names
coef_SL_F3_tab = coef_SL_F3_tab[,-1]
knitr::kable(round(coef_SL_F3_tab, 3))| SL-12 | SL-12 (log) | SL-12 (all) | SL-90 | SL-90 (log) | SL-90 (all) | |
|---|---|---|---|---|---|---|
| bayesglm | 0.063 | 0.102 | 0.104 | 0.039 | 0.003 | 0.042 |
| earth | 0.022 | 0.007 | 0.086 | 0.103 | 0.068 | 0.053 |
| gam | 0.144 | 0.050 | 0.070 | 0.036 | 0.003 | 0.033 |
| gbm | 0.281 | 0.131 | 0.134 | 0.282 | 0.294 | 0.373 |
| glm | 0.072 | 0.109 | 0.093 | 0.039 | 0.003 | 0.020 |
| glmnet | 0.000 | 0.079 | 0.155 | 0.088 | 0.020 | 0.050 |
| ipredbagg | 0.070 | 0.002 | 0.075 | 0.000 | 0.035 | 0.027 |
| nnet | 0.034 | 0.116 | 0.044 | 0.179 | 0.334 | 0.108 |
| polymars | 0.079 | 0.200 | 0.101 | 0.136 | 0.056 | 0.089 |
| randomForest | 0.235 | 0.179 | 0.139 | 0.092 | 0.150 | 0.147 |
| rpart | 0.000 | 0.025 | 0.000 | 0.005 | 0.033 | 0.058 |
| svm | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 |
# Use more descriptive model names
# F2 outcome
coef_SL_F2_tab = merge(sl_dict, coef_SL_F2_tab,
by.x = "short", by.y = 0, sort = FALSE)
rownames(coef_SL_F2_tab) = paste(coef_SL_F2_tab$short, coef_SL_F2_tab$long)
coef_SL_F2_tab$short = NULL; coef_SL_F2_tab$long = NULL
write.csv(round(coef_SL_F2_tab, 3), file = "figs_and_tabs/coef_SL_F2_tab.csv")
# F3 outcome
coef_SL_F3_tab = merge(sl_dict, coef_SL_F3_tab,
by.x = "short", by.y = 0, sort = FALSE)
rownames(coef_SL_F3_tab) = paste(coef_SL_F3_tab$short, coef_SL_F3_tab$long)
coef_SL_F3_tab$short = NULL; coef_SL_F3_tab$long = NULL
write.csv(round(coef_SL_F3_tab, 3), file = "figs_and_tabs/coef_SL_F3_tab.csv")Validate superlearners
We predict the scores for each superlearner model and test dataset. We also calculate the six other fibrosis scores under consideration (APRI, BARD, FIB4, Forns, NFS, and SAFE). The SAFE logistic regression model can be re-fit using the NASH data to obtain SAFE predictions on the probability scale.
# Replicate SAFE score
fit_SAFE = glm(sigfibro2 ~ AGE + BMI40 + DIAB2 + LN_AST +
LN_ALT + LN_GLOB + LN_PLAT,
data = nash, family = binomial)
# Function to calculate APRI, BARD, FIB-4, Forns, NFS, and SAFE scores
calc_other_scores = function(dat, AST_upper = 40) {
# Initialize data frame for storing scores
scores = data.frame(matrix(, nrow=nrow(dat), ncol = 0))
# APRI
scores$APRI = 100 * (dat$AST / AST_upper) / dat$PLAT
# BARD
scores$BARD = as.numeric(
2*((dat$AST / dat$ALT) >= 0.8) + dat$DIAB2 + (dat$BMI >= 28)
)
# FIB-4
scores$FIB4 = (dat$AGE * dat$AST) / (dat$PLAT * sqrt(dat$ALT))
# Forns
scores$Forns = 7.811 - 3.131 * log(dat$PLAT) + 0.781 * log(dat$GGT) +
3.467 * log(dat$AGE) - 0.014 * dat$CHOL
# NFS
scores$NFS = (-1.675) + (0.037 * dat$AGE) + (0.094 * dat$BMI) +
(1.13 * dat$DIAB2) + (0.99 * (dat$AST / dat$ALT)) -
(0.013 * dat$PLAT) - (0.66 * dat$ALB)
# SAFE
scores$SAFE = 2.97 * dat$AGE +
5.99 * dat$BMI40 +
62.85 * dat$DIAB2 +
154.85 * log(dat$AST) -
58.23 * log(dat$ALT) +
195.48 * log(dat$GLOB)-
141.61 * log(dat$PLAT)-
75
# SAFE (probability)
scores$SAFE_prob = predict(fit_SAFE, dat, type = "response")
return(scores)
}
# FLINT scores
scores_flint = data.frame(
# Superlearner with F2 outcome
SL_12_F2 = predict(SL_12_F2, flintX)$pred,
SL_12_log_F2 = predict(SL_12_log_F2, flintX_log)$pred,
SL_12_all_F2 = predict(SL_12_all_F2, flintX_all)$pred,
SL_90_F2 = predict(SL_90_F2, flintX)$pred,
SL_90_log_F2 = predict(SL_90_log_F2, flintX_log)$pred,
SL_90_all_F2 = predict(SL_90_all_F2, flintX_all)$pred,
# Superlearner with F3 outcome
SL_12_F3 = predict(SL_12_F3, flintX)$pred,
SL_12_log_F3 = predict(SL_12_log_F3, flintX_log)$pred,
SL_12_all_F3 = predict(SL_12_all_F3, flintX_all)$pred,
SL_90_F3 = predict(SL_90_F3, flintX)$pred,
SL_90_log_F3 = predict(SL_90_log_F3, flintX_log)$pred,
SL_90_all_F3 = predict(SL_90_all_F3, flintX_all)$pred,
# Other scores
calc_other_scores(flint))
# NHANES-NAFLD scores
scores_nhanes = data.frame(
# Superlearner with F2 outcome
SL_12_F2 = predict(SL_12_F2, nhanesX)$pred,
SL_12_log_F2 = predict(SL_12_log_F2, nhanesX_log)$pred,
SL_12_all_F2 = predict(SL_12_all_F2, nhanesX_all)$pred,
SL_90_F2 = predict(SL_90_F2, nhanesX)$pred,
SL_90_log_F2 = predict(SL_90_log_F2, nhanesX_log)$pred,
SL_90_all_F2 = predict(SL_90_all_F2, nhanesX_all)$pred,
# Superlearner with F3 outcome
SL_12_F3 = predict(SL_12_F3, nhanesX)$pred,
SL_12_log_F3 = predict(SL_12_log_F3, nhanesX_log)$pred,
SL_12_all_F3 = predict(SL_12_all_F3, nhanesX_all)$pred,
SL_90_F3 = predict(SL_90_F3, nhanesX)$pred,
SL_90_log_F3 = predict(SL_90_log_F3, nhanesX_log)$pred,
SL_90_all_F3 = predict(SL_90_all_F3, nhanesX_all)$pred,
# Other scores
calc_other_scores(nhanes))AUCs and ROCs
In each outcome and validation set, we calculate the AUC for each of the predicted superlearner scores as well as the other fibrosis scores. The NHANES-NAFLD AUCs incorporate sampling weights.
# Function to calculate AUC for multiple scores
calc_AUC = function(preds, response) {
sapply(preds, function(x) {
AUC = cvAUC::AUC(x, response)
})
}
# Function to calculate AUC for multiple scores,
# while accounting for sampling weights
calc_AUC_weighted = function(preds, response, weights) {
sapply(preds, function(x) {
AUC = WeightedAUC(WeightedROC(x, response, weights))
})
}
# Column indices for each type of score
SL_F2_idx = 1:6
SL_F3_idx = 7:12
other_idx = 13:18
# F2 outcome
# FLINT
AUC_flint_F2 = calc_AUC(scores_flint[,c(SL_F2_idx, other_idx)], flintY2)
# NHANES-NAFLD
AUC_nhanes_F2 = calc_AUC_weighted(scores_nhanes[,c(SL_F2_idx, other_idx)],
nhanesY2, nhanes_weights)
# F3 outcome
# FLINT
AUC_flint_F3 = calc_AUC(scores_flint[,c(SL_F3_idx, other_idx)], flintY3)
# NHANES-NAFLD
AUC_nhanes_F3 = calc_AUC_weighted(scores_nhanes[,c(SL_F3_idx, other_idx)],
nhanesY3, nhanes_weights)Confidence intervals for the AUCs are obtained by bootstrapping the data 1000 times and taking the 95% percentiles of the bootstrapped AUCs.
# Function to calculate bootstrap AUCs for multiple scores
boot_AUC = function(preds, response, B = 1000) {
run_boot = function(){
# Sample indices with replacement
boot_idx = sample(length(response), replace = TRUE)
# Calculate bootstrap AUC
return(sapply(preds, function(x) {
cvAUC::AUC(x[boot_idx], response[boot_idx])
}))
}
return(t(replicate(B, run_boot())))
}
# Function to calculate bootstrap weighted AUCs for multiple scores
boot_AUC_weighted = function(preds, response, weights, B = 1000) {
run_boot = function(){
# Sample indices with replacement
boot_idx = sample(length(response), replace = TRUE)
# Calculate bootstrap AUC
return(sapply(preds, function(x) {
WeightedAUC(WeightedROC(x[boot_idx], response[boot_idx],
weights[boot_idx]))
}))
}
return(t(replicate(B, run_boot())))
}
# F2 outcome
# FLINT
set.seed(1)
boot_AUC_flint_F2 = boot_AUC(scores_flint[,c(SL_F2_idx, other_idx)], flintY2)
# NHANES-NAFLD
set.seed(1)
boot_AUC_nhanes_F2 = boot_AUC_weighted(scores_nhanes[,c(SL_F2_idx, other_idx)],
nhanesY2, nhanes_weights)
save(boot_AUC_flint_F2, boot_AUC_nhanes_F2,
file = "model_output/boot_AUC_F2.rData")
# F3 outcome
# FLINT
set.seed(1)
boot_AUC_flint_F3 = boot_AUC(scores_flint[,c(SL_F3_idx, other_idx)], flintY3)
# NHANES-NAFLD
set.seed(1)
boot_AUC_nhanes_F3 = boot_AUC_weighted(scores_nhanes[,c(SL_F3_idx, other_idx)],
nhanesY3, nhanes_weights)
save(boot_AUC_flint_F3, boot_AUC_nhanes_F3,
file = "model_output/boot_AUC_F3.rData")We present the AUCs and 95% confidence intervals for each outcome, model, and validation dataset as a table.
# AUCs with 95% CIs
# F2 outcome
AUC_CI_F2_df = cbind(
"FLINT" = AUC_flint_F2,
"FLINT 2.5%" = apply(boot_AUC_flint_F2, 2, quantile, 0.025),
"FLINT 97.5%" = apply(boot_AUC_flint_F2, 2, quantile, 0.975),
"NHANES-NAFLD" = AUC_nhanes_F2,
"NHANES-NAFLD 2.5%" = apply(boot_AUC_nhanes_F2, 2, quantile, 0.025),
"NHANES-NAFLD 97.5%" = apply(boot_AUC_nhanes_F2, 2, quantile, 0.975))
# F3 outcome
AUC_CI_F3_df = cbind(
"FLINT" = AUC_flint_F3,
"FLINT 2.5%" = apply(boot_AUC_flint_F3, 2, quantile, 0.025),
"FLINT 97.5%" = apply(boot_AUC_flint_F3, 2, quantile, 0.975),
"NHANES-NAFLD" = AUC_nhanes_F3,
"NHANES-NAFLD 2.5%" = apply(boot_AUC_nhanes_F3, 2, quantile, 0.025),
"NHANES-NAFLD 97.5%" = apply(boot_AUC_nhanes_F3, 2, quantile, 0.975))
# Formatting as table
AUC_CI_table = data.frame(
matrix(paste0(
round(AUC_CI_F2_df[,seq(1, ncol(AUC_CI_F2_df), by = 3)], 3), " (",
round(AUC_CI_F2_df[,seq(2, ncol(AUC_CI_F2_df), by = 3)], 3), ", ",
round(AUC_CI_F2_df[,seq(3, ncol(AUC_CI_F2_df), by = 3)], 3), ")"),
nrow = nrow(AUC_CI_F2_df), ncol = ncol(AUC_CI_F2_df)/3),
matrix(paste0(
round(AUC_CI_F3_df[,seq(1, ncol(AUC_CI_F3_df), by = 3)], 3), " (",
round(AUC_CI_F3_df[,seq(2, ncol(AUC_CI_F3_df), by = 3)], 3), ", ",
round(AUC_CI_F3_df[,seq(3, ncol(AUC_CI_F3_df), by = 3)], 3), ")"),
nrow = nrow(AUC_CI_F3_df), ncol = ncol(AUC_CI_F3_df)/3)
)
rownames(AUC_CI_table) = c("SL-12", "SL-12 (log)", "SL-12 (all)",
"SL-90", "SL-90 (log)", "SL-90 (all)",
"APRI","BARD", "FIB-4", "Forns", "NFS", "SAFE")
colnames(AUC_CI_table) = paste(c("FLINT", "NHANES-NAFLD"),
rep(c("(F2 or higher)", "(F3 or higher)"),
each = 2))
write.csv(AUC_CI_table, file = "figs_and_tabs/AUC_CI_table.csv")
knitr::kable(AUC_CI_table)| FLINT (F2 or higher) | NHANES-NAFLD (F2 or higher) | FLINT (F3 or higher) | NHANES-NAFLD (F3 or higher) | |
|---|---|---|---|---|
| SL-12 | 0.788 (0.73, 0.836) | 0.711 (0.655, 0.763) | 0.783 (0.72, 0.843) | 0.787 (0.705, 0.852) |
| SL-12 (log) | 0.787 (0.726, 0.833) | 0.738 (0.678, 0.792) | 0.788 (0.726, 0.845) | 0.777 (0.689, 0.846) |
| SL-12 (all) | 0.793 (0.733, 0.841) | 0.71 (0.65, 0.767) | 0.791 (0.729, 0.851) | 0.78 (0.694, 0.855) |
| SL-90 | 0.793 (0.733, 0.84) | 0.719 (0.659, 0.773) | 0.789 (0.725, 0.845) | 0.784 (0.696, 0.852) |
| SL-90 (log) | 0.781 (0.721, 0.829) | 0.705 (0.644, 0.763) | 0.789 (0.724, 0.846) | 0.754 (0.665, 0.836) |
| SL-90 (all) | 0.798 (0.74, 0.846) | 0.686 (0.625, 0.749) | 0.786 (0.724, 0.846) | 0.774 (0.69, 0.846) |
| APRI | 0.716 (0.657, 0.774) | 0.623 (0.56, 0.682) | 0.69 (0.621, 0.759) | 0.669 (0.574, 0.764) |
| BARD | 0.659 (0.595, 0.721) | 0.568 (0.502, 0.637) | 0.656 (0.584, 0.723) | 0.657 (0.53, 0.767) |
| FIB-4 | 0.736 (0.675, 0.792) | 0.59 (0.523, 0.651) | 0.712 (0.641, 0.778) | 0.632 (0.535, 0.726) |
| Forns | 0.65 (0.583, 0.711) | 0.619 (0.559, 0.674) | 0.668 (0.598, 0.732) | 0.634 (0.534, 0.724) |
| NFS | 0.706 (0.644, 0.762) | 0.724 (0.668, 0.782) | 0.683 (0.614, 0.749) | 0.775 (0.696, 0.847) |
| SAFE | 0.795 (0.739, 0.841) | 0.735 (0.674, 0.789) | 0.761 (0.698, 0.818) | 0.803 (0.711, 0.875) |
And as a plot. For brevity, we only show the results from the superlearner with 12 base models and all predictors for the remainder of the analysis.
# Models to present in main text
mod_names_F2 = c("Superlearner-F2", "APRI","BARD", "FIB-4", "Forns", "NFS", "SAFE")
mod_sub_names_F2 = c("SL_12_all_F2",
"APRI", "BARD", "FIB4", "Forns", "NFS", "SAFE")
mod_names_F3 = c("Superlearner-F3", "APRI","BARD", "FIB-4", "Forns", "NFS", "SAFE")
mod_sub_names_F3 = c("SL_12_all_F3",
"APRI", "BARD", "FIB4", "Forns", "NFS", "SAFE")
# Plot AUCs with CIs
AUC_CI_df = rbind(AUC_CI_F2_df[mod_sub_names_F2,],
AUC_CI_F3_df[mod_sub_names_F3,])
rownames(AUC_CI_df) = c(mod_names_F2, mod_names_F3)
p1 = melt(AUC_CI_df[,seq(1, ncol(AUC_CI_df), by = 3)]) %>%
rename(Model = Var1, Data = Var2, est = value) %>%
mutate(lo = melt(AUC_CI_df[,seq(2, ncol(AUC_CI_df), by = 3)])$value,
hi = melt(AUC_CI_df[,seq(3, ncol(AUC_CI_df), by = 3)])$value,
x = rep(1:2, each = 14) + rep(seq(-0.4, 0.4, length = 7), 2) +
rep(c(-0.025, 0.025), each = 7),
Outcome = rep(rep(c("F2", "F3"), each = 7), 2),
Model = factor(Model, levels = c(mod_names_F2[1], mod_names_F3))) %>%
ggplot(aes(x = x, y = est, col = Model, shape = Outcome)) +
geom_pointrange(aes(ymin = lo, ymax = hi), size = 0.5, stroke = 0.5) +
scale_shape_manual(name = "Significant\nfibrosis (data)", values = c(19, 2),
labels = c("F2 or higher", "F3 or higher")) +
xlab("") + ylab("AUC") +
scale_x_continuous(breaks = 1:2,
labels = c("FLINT",
"NHANES-NAFLD")) +
theme(axis.text.x = element_text(size = 10))
ggsave("figs_and_tabs/AUC_CI.png", p1, width = 6.5, height = 4)
p1
# Plot AUCs with CIs for F2 only
p1 = melt(AUC_CI_F2_df[mod_sub_names_F2,seq(1, ncol(AUC_CI_F2_df), by = 3)]) %>%
rename(Model = Var1, Data = Var2, est = value) %>%
mutate(lo = melt(AUC_CI_F2_df[mod_sub_names_F2,seq(2, ncol(AUC_CI_F2_df), by = 3)])$value,
hi = melt(AUC_CI_F2_df[mod_sub_names_F2,seq(3, ncol(AUC_CI_F2_df), by = 3)])$value,
x = rep(1:2, each = 7) + seq(-0.4, 0.4, length = 7),
Model = factor(Model, labels = c("Superlearner", mod_names_F2[-1]))) %>%
ggplot(aes(x = x, y = est, color = Model, shape = Model)) +
geom_pointrange(aes(ymin = lo, ymax = hi), size = 0.5, stroke = 0.5) +
xlab("") + ylab("AUC") +
scale_color_grey() +
scale_shape_manual(values = c(16, 4:6, 0, 8, 17)) +
scale_x_continuous(breaks = 1:2,
labels = c("FLINT",
"NHANES-NAFLD")) +
theme(axis.text.x = element_text(size = 10))
ggsave("figs_and_tabs/AUC_CI_F2.png", p1, width = 6.5, height = 4)
p1
We can also use ROC curves to visualize the each scoring approach ability to discriminate significant fibrosis.
# F2 outcome
# FLINT
p1 = rbind(
do.call(data.frame,
c(roc(flintY2, scores_flint$SL_12_all_F2)[c("sensitivities", "specificities")],
Model = "Superlearner-F2")),
data.frame(sensitivities = NA, specificities = NA, Model = "Superlearner-F3"),
do.call(data.frame,
c(roc(flintY2, scores_flint$APRI)[c("sensitivities", "specificities")],
Model = "APRI")),
do.call(data.frame,
c(roc(flintY2, scores_flint$BARD)[c("sensitivities", "specificities")],
Model = "BARD")),
do.call(data.frame,
c(roc(flintY2, scores_flint$FIB4)[c("sensitivities", "specificities")],
Model = "FIB-4")),
do.call(data.frame,
c(roc(flintY2, scores_flint$Forns)[c("sensitivities", "specificities")],
Model = "Forns")),
do.call(data.frame,
c(roc(flintY2, scores_flint$NFS)[c("sensitivities", "specificities")],
Model = "NFS")),
do.call(data.frame,
c(roc(flintY2, scores_flint$SAFE)[c("sensitivities", "specificities")],
Model = "SAFE"))) %>%
mutate(Model = fct_relevel(Model, c(mod_names_F2[1], mod_names_F3))) %>%
ggplot(aes(x = specificities, y = sensitivities, col = Model)) +
geom_path() +
scale_x_reverse() +
scale_color_discrete(breaks = mod_names_F2) +
geom_abline(intercept = 1, slope = 1, col = "grey") +
xlab("Specificity") + ylab("Sensitivity") +
ggtitle("FLINT")
# NHANES-NAFLD
p2 = rbind(
cbind(WeightedROC(scores_nhanes$SL_12_all_F2, nhanesY2, nhanes_weights),
Model = "Superlearner-F2"),
cbind(data.frame(TPR = NA, FPR = NA, threshold = NA, FN = NA, FP = NA),
Model = "Superlearner-F3"),
cbind(WeightedROC(scores_nhanes$APRI, nhanesY2, nhanes_weights),
Model = "APRI"),
cbind(WeightedROC(scores_nhanes$BARD, nhanesY2, nhanes_weights),
Model = "BARD"),
cbind(WeightedROC(scores_nhanes$FIB4, nhanesY2, nhanes_weights),
Model = "FIB-4"),
cbind(WeightedROC(scores_nhanes$Forns, nhanesY2, nhanes_weights),
Model = "Forns"),
cbind(WeightedROC(scores_nhanes$NFS, nhanesY2, nhanes_weights),
Model = "NFS"),
cbind(WeightedROC(scores_nhanes$SAFE, nhanesY2, nhanes_weights),
Model = "SAFE")) %>%
mutate(Model = fct_relevel(Model, c(mod_names_F2[1], mod_names_F3))) %>%
ggplot(aes(x = 1 - FPR, y = TPR, col = Model)) +
geom_path() +
scale_x_reverse() +
scale_color_discrete(breaks = mod_names_F2) +
geom_abline(intercept = 1, slope = 1, col = "grey") +
xlab("Specificity") + ylab("Sensitivity") +
ggtitle("NHANES-NAFLD")
# Legend
p0 = gtable_filter(ggplot_gtable(ggplot_build(
p2 + theme(legend.margin=margin(l = 0, r = 2, unit='cm')))), "guide-box")
# Put all ROC plots and legend together
png("figs_and_tabs/ROC_F2.png", width = 1100, height = 500, res = 140)
plot_grid(p1 + guides(col = "none"),
p2 + guides(col = "none"),
p0, ncol = 3, rel_widths = c(3/7, 3/7, 1/7))
invisible(dev.off())
plot_grid(p1 + guides(col = "none"),
p2 + guides(col = "none"),
p0, ncol = 3, rel_widths = c(3/7, 3/7, 1/7))
# F2 outcome, black and white
# FLINT
p1 = rbind(
do.call(data.frame,
c(roc(flintY2, scores_flint$SL_12_all_F2)[c("sensitivities", "specificities")],
Model = "Superlearner")),
do.call(data.frame,
c(roc(flintY2, scores_flint$APRI)[c("sensitivities", "specificities")],
Model = "APRI")),
do.call(data.frame,
c(roc(flintY2, scores_flint$BARD)[c("sensitivities", "specificities")],
Model = "BARD")),
do.call(data.frame,
c(roc(flintY2, scores_flint$FIB4)[c("sensitivities", "specificities")],
Model = "FIB-4")),
do.call(data.frame,
c(roc(flintY2, scores_flint$Forns)[c("sensitivities", "specificities")],
Model = "Forns")),
do.call(data.frame,
c(roc(flintY2, scores_flint$NFS)[c("sensitivities", "specificities")],
Model = "NFS")),
do.call(data.frame,
c(roc(flintY2, scores_flint$SAFE)[c("sensitivities", "specificities")],
Model = "SAFE"))) %>%
mutate(Model = factor(Model, c("Superlearner", mod_names_F2[-1]))) %>%
ggplot(aes(x = specificities, y = sensitivities, col = Model, lty = Model)) +
geom_path() +
scale_color_manual(name = "Model",
values = c(rep(c("black", "grey50"), each = 3), "grey80"),
labels = c("Superlearner", mod_names_F2[-1])) +
scale_linetype_manual(name = "Model",
values = rep(c("solid", "dashed", "dotted"), length = 7),
labels = c("Superlearner", mod_names_F2[-1])) +
scale_x_reverse() +
geom_abline(intercept = 1, slope = 1, col = "grey") +
xlab("Specificity") + ylab("Sensitivity") +
ggtitle("FLINT")
# NHANES-NAFLD
p2 = rbind(
cbind(WeightedROC(scores_nhanes$SL_12_all_F2, nhanesY2, nhanes_weights),
Model = "Superlearner"),
cbind(WeightedROC(scores_nhanes$APRI, nhanesY2, nhanes_weights),
Model = "APRI"),
cbind(WeightedROC(scores_nhanes$BARD, nhanesY2, nhanes_weights),
Model = "BARD"),
cbind(WeightedROC(scores_nhanes$FIB4, nhanesY2, nhanes_weights),
Model = "FIB-4"),
cbind(WeightedROC(scores_nhanes$Forns, nhanesY2, nhanes_weights),
Model = "Forns"),
cbind(WeightedROC(scores_nhanes$NFS, nhanesY2, nhanes_weights),
Model = "NFS"),
cbind(WeightedROC(scores_nhanes$SAFE, nhanesY2, nhanes_weights),
Model = "SAFE")) %>%
mutate(Model = factor(Model, c("Superlearner", mod_names_F2[-1]))) %>%
ggplot(aes(x = 1 - FPR, y = TPR, col = Model, lty = Model)) +
geom_path() +
scale_color_manual(name = "Model",
values = c(rep(c("black", "grey50"), each = 3), "grey80"),
labels = c("Superlearner", mod_names_F2[-1])) +
scale_linetype_manual(name = "Model",
values = rep(c("solid", "dashed", "dotted"), length = 7),
labels = c("Superlearner", mod_names_F2[-1])) +
scale_x_reverse() +
geom_abline(intercept = 1, slope = 1, col = "grey") +
xlab("Specificity") + ylab("Sensitivity") +
ggtitle("NHANES-NAFLD")
# Legend
p0 = gtable_filter(ggplot_gtable(ggplot_build(
p2 + theme(legend.margin=margin(l = 0, r = 2, unit='cm')))), "guide-box")
# Put all ROC plots and legend together
png("figs_and_tabs/ROC_F2_bw.png", width = 1100, height = 500, res = 140)
plot_grid(p1 + guides(col = "none", lty = "none"),
p2 + guides(col = "none", lty = "none"),
p0, ncol = 3, rel_widths = c(3/7, 3/7, 1/7))
invisible(dev.off())
plot_grid(p1 + guides(col = "none", lty = "none"),
p2 + guides(col = "none", lty = "none"),
p0, ncol = 3, rel_widths = c(3/7, 3/7, 1/7))
# F3 outcome
# FLINT
p1 = rbind(
data.frame(sensitivities = NA, specificities = NA, Model = "Superlearner-F2"),
do.call(data.frame,
c(roc(flintY3, scores_flint$SL_12_all_F3)[c("sensitivities", "specificities")],
Model = "Superlearner-F3")),
do.call(data.frame,
c(roc(flintY3, scores_flint$APRI)[c("sensitivities", "specificities")],
Model = "APRI")),
do.call(data.frame,
c(roc(flintY3, scores_flint$BARD)[c("sensitivities", "specificities")],
Model = "BARD")),
do.call(data.frame,
c(roc(flintY3, scores_flint$FIB4)[c("sensitivities", "specificities")],
Model = "FIB-4")),
do.call(data.frame,
c(roc(flintY3, scores_flint$Forns)[c("sensitivities", "specificities")],
Model = "Forns")),
do.call(data.frame,
c(roc(flintY3, scores_flint$NFS)[c("sensitivities", "specificities")],
Model = "NFS")),
do.call(data.frame,
c(roc(flintY3, scores_flint$SAFE)[c("sensitivities", "specificities")],
Model = "SAFE"))) %>%
mutate(Model = fct_relevel(Model, c(mod_names_F2[1], mod_names_F3))) %>%
ggplot(aes(x = specificities, y = sensitivities, col = Model)) +
geom_path() +
scale_x_reverse() +
scale_color_discrete(breaks = mod_names_F3) +
geom_abline(intercept = 1, slope = 1, col = "grey") +
xlab("Specificity") + ylab("Sensitivity") +
ggtitle("FLINT")
# NHANES-NAFLD
p2 = rbind(
cbind(data.frame(TPR = NA, FPR = NA, threshold = NA, FN = NA, FP = NA),
Model = "Superlearner-F2"),
cbind(WeightedROC(scores_nhanes$SL_12_all_F3, nhanesY3, nhanes_weights),
Model = "Superlearner-F3"),
cbind(WeightedROC(scores_nhanes$APRI, nhanesY3, nhanes_weights),
Model = "APRI"),
cbind(WeightedROC(scores_nhanes$BARD, nhanesY3, nhanes_weights),
Model = "BARD"),
cbind(WeightedROC(scores_nhanes$FIB4, nhanesY3, nhanes_weights),
Model = "FIB-4"),
cbind(WeightedROC(scores_nhanes$Forns, nhanesY3, nhanes_weights),
Model = "Forns"),
cbind(WeightedROC(scores_nhanes$NFS, nhanesY3, nhanes_weights),
Model = "NFS"),
cbind(WeightedROC(scores_nhanes$SAFE, nhanesY3, nhanes_weights),
Model = "SAFE")) %>%
mutate(Model = fct_relevel(Model, c(mod_names_F2[1], mod_names_F3))) %>%
ggplot(aes(x = 1 - FPR, y = TPR, col = Model)) +
geom_path() +
scale_x_reverse() +
scale_color_discrete(breaks = mod_names_F3) +
geom_abline(intercept = 1, slope = 1, col = "grey") +
xlab("Specificity") + ylab("Sensitivity") +
ggtitle("NHANES-NAFLD")
# Legend
p0 = gtable_filter(ggplot_gtable(ggplot_build(
p2 + theme(legend.margin=margin(l = 0, r = 2, unit='cm')))), "guide-box")
# Put all ROC plots and legend together
png("figs_and_tabs/ROC_F3.png", width = 1100, height = 500, res = 140)
plot_grid(p1 + guides(col = "none"),
p2 + guides(col = "none"),
p0, ncol = 3, rel_widths = c(3/7, 3/7, 1/7))
invisible(dev.off())
plot_grid(p1 + guides(col = "none"),
p2 + guides(col = "none"),
p0, ncol = 3, rel_widths = c(3/7, 3/7, 1/7))
Calibration and Brier scores
Plots of the observed risk vs. predicted probabilities for each of the superlearners and the refitted SAFE model, with loess smoothers.
# Column indices for scores that are probabilities
prob_F2_idx = c(SL_F2_idx, 19)
prob_F3_idx = c(SL_F3_idx, 19)
# FLINT
# F2 or higher
data.frame(melt(scores_flint[,prob_F2_idx]), Observed = flintY2) %>%
rename(Model = variable, `Predicted probability` = value) %>%
mutate(Model = factor(Model,
labels = c("SL-12", "SL-12 (log)", "SL-12 (all)",
"SL-90", "SL-90 (log)", "SL-90 (all)", "SAFE"))) %>%
ggplot(aes(x = `Predicted probability`, y = Observed, color = Model)) +
geom_smooth(method = "loess", se = FALSE, lwd = 1) +
geom_abline(intercept = 0, slope = 1, color = "grey") +
ggtitle("FLINT (F2 or higher)")
# F3 or higher
data.frame(melt(scores_flint[,prob_F3_idx]), Observed = flintY3) %>%
rename(Model = variable, `Predicted probability` = value) %>%
mutate(Model = factor(Model,
labels = c("SL-12", "SL-12 (log)", "SL-12 (all)",
"SL-90", "SL-90 (log)", "SL-90 (all)", "SAFE"))) %>%
ggplot(aes(x = `Predicted probability`, y = Observed, color = Model)) +
geom_smooth(method = "loess", se = FALSE, lwd = 1) +
geom_abline(intercept = 0, slope = 1, color = "grey") +
ggtitle("FLINT (F3 or higher)")
# NHANES-NAFLD
# F2 or higher
data.frame(melt(scores_nhanes[,prob_F2_idx]), Observed = nhanesY2) %>%
rename(Model = variable, `Predicted probability` = value) %>%
mutate(Model = factor(Model,
labels = c("SL-12", "SL-12 (log)", "SL-12 (all)",
"SL-90", "SL-90 (log)", "SL-90 (all)", "SAFE"))) %>%
ggplot(aes(x = `Predicted probability`, y = Observed, color = Model)) +
geom_smooth(method = "loess", se = FALSE, lwd = 1) +
geom_abline(intercept = 0, slope = 1, color = "grey") +
ggtitle("NHANES-NAFLD (F2 or higher)")
# F3 or higher
data.frame(melt(scores_nhanes[,prob_F3_idx]), Observed = nhanesY3) %>%
rename(Model = variable, `Predicted probability` = value) %>%
mutate(Model = factor(Model,
labels = c("SL-12", "SL-12 (log)", "SL-12 (all)",
"SL-90", "SL-90 (log)", "SL-90 (all)", "SAFE"))) %>%
ggplot(aes(x = `Predicted probability`, y = Observed, color = Model)) +
geom_smooth(method = "loess", se = FALSE, lwd = 1) +
geom_abline(intercept = 0, slope = 1, color = "grey") +
ggtitle("NHANES-NAFLD (F3 or higher)")
Again, but restricting to the SL-12 (all) and SAFE models.
# Column indices for scores that are probabilities
prob_F2_idx = c(SL_F2_idx, 19)
prob_F3_idx = c(SL_F3_idx, 19)
# FLINT
# F2 or higher
data.frame(melt(scores_flint[,c("SL_90_all_F2", "SAFE_prob")]), Observed = flintY2) %>%
rename(Model = variable, `Predicted probability` = value) %>%
mutate(Model = factor(Model,
labels = c("Superlearner", "SAFE"))) %>%
ggplot(aes(x = `Predicted probability`, y = Observed, color = Model)) +
geom_smooth(method = "loess", se = FALSE, lwd = 1) +
geom_abline(intercept = 0, slope = 1, color = "grey") +
ggtitle("FLINT (F2 or higher)")
# F3 or higher
data.frame(melt(scores_flint[,c("SL_90_all_F3", "SAFE_prob")]), Observed = flintY3) %>%
rename(Model = variable, `Predicted probability` = value) %>%
mutate(Model = factor(Model,
labels = c("Superlearner", "SAFE"))) %>%
ggplot(aes(x = `Predicted probability`, y = Observed, color = Model)) +
geom_smooth(method = "loess", se = FALSE, lwd = 1) +
geom_abline(intercept = 0, slope = 1, color = "grey") +
ggtitle("FLINT (F3 or higher)")
# NHANES-NAFLD
# F2 or higher
data.frame(melt(scores_nhanes[,c("SL_90_all_F2", "SAFE_prob")]), Observed = nhanesY2) %>%
rename(Model = variable, `Predicted probability` = value) %>%
mutate(Model = factor(Model,
labels = c("Superlearner", "SAFE"))) %>%
ggplot(aes(x = `Predicted probability`, y = Observed, color = Model)) +
geom_smooth(method = "loess", se = FALSE, lwd = 1) +
geom_abline(intercept = 0, slope = 1, color = "grey") +
ggtitle("NHANES-NAFLD (F2 or higher)")
# F3 or higher
data.frame(melt(scores_nhanes[,c("SL_90_all_F3", "SAFE_prob")]), Observed = nhanesY3) %>%
rename(Model = variable, `Predicted probability` = value) %>%
mutate(Model = factor(Model,
labels = c("Superlearner", "SAFE"))) %>%
ggplot(aes(x = `Predicted probability`, y = Observed, color = Model)) +
geom_smooth(method = "loess", se = FALSE, lwd = 1) +
geom_abline(intercept = 0, slope = 1, color = "grey") +
ggtitle("NHANES-NAFLD (F3 or higher)")
In each outcome and validation set, we calculate the overall calibration (observed / expected number of patients with significant fibrosis) and Brier score for each of the predicted superlearner scores and the SAFE probabilities. Using the SAFE probabilities, we group the patients in each data set into deciles and calculate the calibration within each decile. SAFE The NHANES-NAFLD performance metrics incorporate sampling weights.
# Function to calculate calibration (observed / expected) for multiple scores
calc_OE = function(preds, response) {
sapply(preds, function(x) {
OE = sum(response) / sum(x)
})
}
# Function to calculate calibration (observed / expected) for multiple scores,
# while accounting for sampling weights
calc_OE_weighted = function(preds, response, weights) {
sapply(preds, function(x) {
OE = weighted.mean(response, weights) / weighted.mean(x, weights)
})
}
# Function to calculate Brier score for multiple scores
calc_Brier = function(preds, response) {
sapply(preds, function(x) {
Brier = mean((x - response)^2)
})
}
# Function to calculate Brier score for multiple scores,
# while accounting for sampling weights
calc_Brier_weighted = function(preds, response, weights) {
sapply(preds, function(x) {
Brier = weighted.mean((x - response)^2, weights)
})
}
# SAFE deciles for each validation set
# FLINT
SAFE_dec_flint_thresh = quantile(scores_flint$SAFE_prob, seq(0.1, 0.9, by = 0.1))
SAFE_dec_flint =
as.numeric(cut(scores_flint$SAFE_prob,
breaks = c(-Inf, as.numeric(SAFE_dec_flint_thresh), Inf),
labels = c(1:(length(SAFE_dec_flint_thresh)+1))))
# NHANES
SAFE_dec_nhanes_thresh = quantile(scores_nhanes$SAFE_prob, seq(0.1, 0.9, by = 0.1))
SAFE_dec_nhanes =
as.numeric(cut(scores_nhanes$SAFE_prob,
breaks = c(-Inf, as.numeric(SAFE_dec_nhanes_thresh), Inf),
labels = c(1:(length(SAFE_dec_nhanes_thresh)+1))))
# Calibration
# F2 outcome
# FLINT
OE_flint_F2 = calc_OE(scores_flint[,prob_F2_idx], flintY2)
# NHANES-NAFLD
OE_nhanes_F2 = calc_OE_weighted(scores_nhanes[,prob_F2_idx],
nhanesY2, nhanes_weights)
# F3 outcome
# FLINT
OE_flint_F3 = calc_OE(scores_flint[,prob_F3_idx], flintY3)
# NHANES-NAFLD
OE_nhanes_F3 = calc_OE_weighted(scores_nhanes[,prob_F3_idx],
nhanesY3, nhanes_weights)
# Calibration by SAFE decile
# F2 outcome
# FLINT
OE_decile_flint_F2 = t(sapply(1:10, function(i) {
calc_OE(scores_flint[SAFE_dec_flint == i, prob_F2_idx],
flintY2[SAFE_dec_flint==i])
}))
# NHANES-NAFLD
OE_decile_nhanes_F2 = t(sapply(1:10, function(i) {
calc_OE_weighted(scores_nhanes[SAFE_dec_nhanes == i, prob_F2_idx],
nhanesY2[SAFE_dec_nhanes == i],
nhanes_weights[SAFE_dec_nhanes == i])
}))
# F3 outcome
# FLINT
OE_decile_flint_F3 = t(sapply(1:10, function(i) {
calc_OE(scores_flint[SAFE_dec_flint == i, prob_F3_idx],
flintY3[SAFE_dec_flint==i])
}))
# NHANES-NAFLD
OE_decile_nhanes_F3 = t(sapply(1:10, function(i) {
calc_OE_weighted(scores_nhanes[SAFE_dec_nhanes == i, prob_F3_idx],
nhanesY3[SAFE_dec_nhanes == i],
nhanes_weights[SAFE_dec_nhanes == i])
}))
# Brier score
# F2 outcome
# FLINT
Brier_flint_F2 = calc_Brier(scores_flint[,prob_F2_idx], flintY2)
# NHANES-NAFLD
Brier_nhanes_F2 = calc_Brier_weighted(scores_nhanes[,prob_F2_idx],
nhanesY2, nhanes_weights)
# F3 outcome
# FLINT
Brier_flint_F3 = calc_Brier(scores_flint[,prob_F3_idx], flintY3)
# NHANES-NAFLD
Brier_nhanes_F3 = calc_Brier_weighted(scores_nhanes[,prob_F3_idx],
nhanesY3, nhanes_weights)Confidence intervals for calibration and Brier score are obtained by bootstrapping the data 1000 times and taking the 95% percentiles of the bootstrapped metrics.
# Function to calculate bootstrap calibration (observed / expected) for
# multiple scores
boot_OE = function(preds, response, B = 1000) {
run_boot = function(){
# Sample indices with replacement
boot_idx = sample(length(response), replace = TRUE)
# Calculate bootstrap AUC
return(sapply(preds, function(x) {
sum(response[boot_idx]) / sum(x[boot_idx])
}))
}
return(t(replicate(B, run_boot())))
}
# Function to calculate bootstrap weighted calibration (observed / expected) for
# multiple scores
boot_OE_weighted = function(preds, response, weights, B = 1000) {
run_boot = function(){
# Sample indices with replacement
boot_idx = sample(length(response), replace = TRUE)
# Calculate bootstrap AUC
return(sapply(preds, function(x) {
weighted.mean(response[boot_idx], weights[boot_idx]) /
weighted.mean(x[boot_idx], weights[boot_idx])
}))
}
return(t(replicate(B, run_boot())))
}
# Function to calculate bootstrap Brier score for multiple scores
boot_Brier = function(preds, response, B = 1000) {
run_boot = function(){
# Sample indices with replacement
boot_idx = sample(length(response), replace = TRUE)
# Calculate bootstrap AUC
return(sapply(preds, function(x) {
mean((x[boot_idx] - response[boot_idx])^2)
}))
}
return(t(replicate(B, run_boot())))
}
# Function to calculate bootstrap weighted Brier score for multiple scores
boot_Brier_weighted = function(preds, response, weights, B = 1000) {
run_boot = function(){
# Sample indices with replacement
boot_idx = sample(length(response), replace = TRUE)
# Calculate bootstrap AUC
return(sapply(preds, function(x) {
weighted.mean((x[boot_idx] - response[boot_idx])^2, weights[boot_idx])
}))
}
return(t(replicate(B, run_boot())))
}
# Calibration
# F2 outcome
# FLINT
set.seed(1)
boot_OE_flint_F2 = boot_OE(scores_flint[,prob_F2_idx], flintY2)
# NHANES-NAFLD
set.seed(1)
boot_OE_nhanes_F2 = boot_OE_weighted(scores_nhanes[,prob_F2_idx],
nhanesY2, nhanes_weights)
# F3 outcome
# FLINT
set.seed(1)
boot_OE_flint_F3 = boot_OE(scores_flint[,prob_F3_idx], flintY3)
# NHANES-NAFLD
set.seed(1)
boot_OE_nhanes_F3 = boot_OE_weighted(scores_nhanes[,prob_F3_idx],
nhanesY3, nhanes_weights)
# Calibration by SAFE decile
# F2 outcome
# FLINT
set.seed(1)
boot_OE_decile_flint_F2 = lapply(1:10, function(i) {
boot_OE(scores_flint[SAFE_dec_flint == i, prob_F2_idx],
flintY2[SAFE_dec_flint==i])
})
# NHANES-NAFLD
set.seed(1)
boot_OE_decile_nhanes_F2 = lapply(1:10, function(i) {
boot_OE_weighted(scores_nhanes[SAFE_dec_nhanes == i, prob_F2_idx],
nhanesY2[SAFE_dec_nhanes == i],
nhanes_weights[SAFE_dec_nhanes == i])
})
# F3 outcome
# FLINT
set.seed(1)
boot_OE_decile_flint_F3 = lapply(1:10, function(i) {
boot_OE(scores_flint[SAFE_dec_flint == i, prob_F3_idx],
flintY3[SAFE_dec_flint==i])
})
# NHANES-NAFLD
set.seed(1)
boot_OE_decile_nhanes_F3 = lapply(1:10, function(i) {
boot_OE_weighted(scores_nhanes[SAFE_dec_nhanes == i, prob_F3_idx],
nhanesY3[SAFE_dec_nhanes == i],
nhanes_weights[SAFE_dec_nhanes == i])
})
# Brier score
# F2 outcome
# FLINT
set.seed(1)
boot_Brier_flint_F2 = boot_Brier(scores_flint[,prob_F2_idx], flintY2)
# NHANES-NAFLD
set.seed(1)
boot_Brier_nhanes_F2 = boot_Brier_weighted(scores_nhanes[,prob_F2_idx],
nhanesY2, nhanes_weights)
# F3 outcome
# FLINT
set.seed(1)
boot_Brier_flint_F3 = boot_Brier(scores_flint[,prob_F3_idx], flintY3)
# NHANES-NAFLD
set.seed(1)
boot_Brier_nhanes_F3 = boot_Brier_weighted(scores_nhanes[,prob_F3_idx],
nhanesY3, nhanes_weights)
save(boot_OE_flint_F2, boot_OE_nhanes_F2,
boot_OE_decile_flint_F2, boot_OE_decile_nhanes_F2,
boot_Brier_flint_F2, boot_Brier_nhanes_F2,
file = "model_output/boot_OE_Brier_F2.rData")
save(boot_OE_flint_F3, boot_OE_nhanes_F3,
boot_OE_decile_flint_F3, boot_OE_decile_nhanes_F3,
boot_Brier_flint_F3, boot_Brier_nhanes_F3,
file = "model_output/boot_OE_Brier_F3.rData")We present the calibration and Brier scores with 95% confidence intervals for each outcome, model, and validation dataset as a table.
# Calibration (observed / expected) with 95% CIs
# F2 outcome
OE_CI_F2_df = cbind(
"FLINT" = OE_flint_F2,
"FLINT 2.5%" = apply(boot_OE_flint_F2, 2, quantile, 0.025),
"FLINT 97.5%" = apply(boot_OE_flint_F2, 2, quantile, 0.975),
"NHANES-NAFLD" = OE_nhanes_F2,
"NHANES-NAFLD 2.5%" = apply(boot_OE_nhanes_F2, 2, quantile, 0.025),
"NHANES-NAFLD 97.5%" = apply(boot_OE_nhanes_F2, 2, quantile, 0.975))
# F3 outcome
OE_CI_F3_df = cbind(
"FLINT" = OE_flint_F3,
"FLINT 2.5%" = apply(boot_OE_flint_F3, 2, quantile, 0.025),
"FLINT 97.5%" = apply(boot_OE_flint_F3, 2, quantile, 0.975),
"NHANES-NAFLD" = OE_nhanes_F3,
"NHANES-NAFLD 2.5%" = apply(boot_OE_nhanes_F3, 2, quantile, 0.025),
"NHANES-NAFLD 97.5%" = apply(boot_OE_nhanes_F3, 2, quantile, 0.975))
# Formatting as table
OE_CI_table = data.frame(
matrix(paste0(
round(OE_CI_F2_df[,seq(1, ncol(OE_CI_F2_df), by = 3)], 3), " (",
round(OE_CI_F2_df[,seq(2, ncol(OE_CI_F2_df), by = 3)], 3), ", ",
round(OE_CI_F2_df[,seq(3, ncol(OE_CI_F2_df), by = 3)], 3), ")"),
nrow = nrow(OE_CI_F2_df), ncol = ncol(OE_CI_F2_df)/3),
matrix(paste0(
round(OE_CI_F3_df[,seq(1, ncol(OE_CI_F3_df), by = 3)], 3), " (",
round(OE_CI_F3_df[,seq(2, ncol(OE_CI_F3_df), by = 3)], 3), ", ",
round(OE_CI_F3_df[,seq(3, ncol(OE_CI_F3_df), by = 3)], 3), ")"),
nrow = nrow(OE_CI_F3_df), ncol = ncol(OE_CI_F3_df)/3)
)
rownames(OE_CI_table) = c("SL-12", "SL-12 (log)", "SL-12 (all)",
"SL-90", "SL-90 (log)", "SL-90 (all)",
"SAFE")
colnames(OE_CI_table) = paste(c("FLINT", "NHANES-NAFLD"),
rep(c("(F2 or higher)", "(F3 or higher)"),
each = 2))
write.csv(OE_CI_table, file = "figs_and_tabs/OE_CI_table.csv")
knitr::kable(OE_CI_table)| FLINT (F2 or higher) | NHANES-NAFLD (F2 or higher) | FLINT (F3 or higher) | NHANES-NAFLD (F3 or higher) | |
|---|---|---|---|---|
| SL-12 | 1.037 (0.941, 1.12) | 0.413 (0.339, 0.488) | 1.07 (0.909, 1.23) | 0.378 (0.251, 0.516) |
| SL-12 (log) | 1.038 (0.941, 1.121) | 0.461 (0.379, 0.542) | 1.196 (1.019, 1.374) | 0.428 (0.286, 0.583) |
| SL-12 (all) | 1.042 (0.945, 1.124) | 0.426 (0.348, 0.502) | 1.04 (0.882, 1.193) | 0.418 (0.278, 0.568) |
| SL-90 | 1.058 (0.959, 1.141) | 0.427 (0.349, 0.504) | 1.261 (1.072, 1.449) | 0.451 (0.299, 0.615) |
| SL-90 (log) | 1.026 (0.929, 1.107) | 0.437 (0.358, 0.515) | 1.356 (1.15, 1.552) | 0.481 (0.316, 0.658) |
| SL-90 (all) | 1.057 (0.959, 1.141) | 0.407 (0.333, 0.481) | 1.104 (0.939, 1.264) | 0.417 (0.278, 0.571) |
| SAFE | 1.045 (0.949, 1.129) | 0.458 (0.376, 0.54) | 0.591 (0.498, 0.686) | 0.178 (0.119, 0.242) |
# Brier scores with 95% CIs
# F2 outcome
Brier_CI_F2_df = cbind(
"FLINT" = Brier_flint_F2,
"FLINT 2.5%" = apply(boot_Brier_flint_F2, 2, quantile, 0.025),
"FLINT 97.5%" = apply(boot_Brier_flint_F2, 2, quantile, 0.975),
"NHANES-NAFLD" = Brier_nhanes_F2,
"NHANES-NAFLD 2.5%" = apply(boot_Brier_nhanes_F2, 2, quantile, 0.025),
"NHANES-NAFLD 97.5%" = apply(boot_Brier_nhanes_F2, 2, quantile, 0.975))
# F3 outcome
Brier_CI_F3_df = cbind(
"FLINT" = Brier_flint_F3,
"FLINT 2.5%" = apply(boot_Brier_flint_F3, 2, quantile, 0.025),
"FLINT 97.5%" = apply(boot_Brier_flint_F3, 2, quantile, 0.975),
"NHANES-NAFLD" = Brier_nhanes_F3,
"NHANES-NAFLD 2.5%" = apply(boot_Brier_nhanes_F3, 2, quantile, 0.025),
"NHANES-NAFLD 97.5%" = apply(boot_Brier_nhanes_F3, 2, quantile, 0.975))
# Formatting as table
Brier_CI_table = data.frame(
matrix(paste0(
round(Brier_CI_F2_df[,seq(1, ncol(Brier_CI_F2_df), by = 3)], 3), " (",
round(Brier_CI_F2_df[,seq(2, ncol(Brier_CI_F2_df), by = 3)], 3), ", ",
round(Brier_CI_F2_df[,seq(3, ncol(Brier_CI_F2_df), by = 3)], 3), ")"),
nrow = nrow(Brier_CI_F2_df), ncol = ncol(Brier_CI_F2_df)/3),
matrix(paste0(
round(Brier_CI_F3_df[,seq(1, ncol(Brier_CI_F3_df), by = 3)], 3), " (",
round(Brier_CI_F3_df[,seq(2, ncol(Brier_CI_F3_df), by = 3)], 3), ", ",
round(Brier_CI_F3_df[,seq(3, ncol(Brier_CI_F3_df), by = 3)], 3), ")"),
nrow = nrow(Brier_CI_F3_df), ncol = ncol(Brier_CI_F3_df)/3)
)
rownames(Brier_CI_table) = c("SL-12", "SL-12 (log)", "SL-12 (all)",
"SL-90", "SL-90 (log)", "SL-90 (all)",
"SAFE")
colnames(Brier_CI_table) = paste(c("FLINT", "NHANES-NAFLD"),
rep(c("(F2 or higher)", "(F3 or higher)"),
each = 2))
write.csv(Brier_CI_table, file = "figs_and_tabs/Brier_CI_table.csv")
knitr::kable(Brier_CI_table)| FLINT (F2 or higher) | NHANES-NAFLD (F2 or higher) | FLINT (F3 or higher) | NHANES-NAFLD (F3 or higher) | |
|---|---|---|---|---|
| SL-12 | 0.187 (0.171, 0.206) | 0.172 (0.16, 0.185) | 0.171 (0.148, 0.194) | 0.064 (0.052, 0.077) |
| SL-12 (log) | 0.186 (0.171, 0.205) | 0.148 (0.137, 0.159) | 0.175 (0.15, 0.2) | 0.059 (0.046, 0.073) |
| SL-12 (all) | 0.185 (0.169, 0.203) | 0.164 (0.153, 0.176) | 0.169 (0.145, 0.191) | 0.059 (0.046, 0.073) |
| SL-90 | 0.185 (0.17, 0.203) | 0.163 (0.152, 0.175) | 0.177 (0.15, 0.201) | 0.057 (0.044, 0.071) |
| SL-90 (log) | 0.187 (0.17, 0.206) | 0.161 (0.149, 0.173) | 0.181 (0.153, 0.209) | 0.056 (0.041, 0.072) |
| SL-90 (all) | 0.184 (0.169, 0.202) | 0.177 (0.164, 0.19) | 0.171 (0.147, 0.195) | 0.06 (0.047, 0.074) |
| SAFE | 0.182 (0.163, 0.203) | 0.157 (0.144, 0.17) | 0.239 (0.214, 0.264) | 0.144 (0.133, 0.156) |
This table summarizes the calibration by SAFE score decile for each outcome, model, and validation dataset.
# Calibration (observed / expected) with 95% CIs
# F2 outcome
OE_decile_CI_F2_df = cbind(
"FLINT" = as.numeric(OE_decile_flint_F2),
"FLINT 2.5%" = as.numeric(t(sapply(1:10, function(i) {
apply(boot_OE_decile_flint_F2[[i]], 2, quantile, 0.025)
}))),
"FLINT 97.5%" = as.numeric(t(sapply(1:10, function(i) {
apply(boot_OE_decile_flint_F2[[i]], 2, quantile, 0.975)
}))),
"NHANES-NAFLD" = as.numeric(OE_decile_nhanes_F2),
"NHANES-NAFLD 2.5%" = as.numeric(t(sapply(1:10, function(i) {
apply(boot_OE_decile_nhanes_F2[[i]], 2, quantile, 0.025)
}))),
"NHANES-NAFLD 97.5%" = as.numeric(t(sapply(1:10, function(i) {
apply(boot_OE_decile_nhanes_F2[[i]], 2, quantile, 0.975)
}))))
# F3 outcome
OE_decile_CI_F3_df = cbind(
"FLINT" = as.numeric(OE_decile_flint_F3),
"FLINT 2.5%" = as.numeric(t(sapply(1:10, function(i) {
apply(boot_OE_decile_flint_F3[[i]], 2, quantile, 0.025)
}))),
"FLINT 97.5%" = as.numeric(t(sapply(1:10, function(i) {
apply(boot_OE_decile_flint_F3[[i]], 2, quantile, 0.975)
}))),
"NHANES-NAFLD" = as.numeric(OE_decile_nhanes_F3),
"NHANES-NAFLD 2.5%" = as.numeric(t(sapply(1:10, function(i) {
apply(boot_OE_decile_nhanes_F3[[i]], 2, quantile, 0.025)
}))),
"NHANES-NAFLD 97.5%" = as.numeric(t(sapply(1:10, function(i) {
apply(boot_OE_decile_nhanes_F3[[i]], 2, quantile, 0.975)
}))))
# Formatting as table
OE_decile_CI_table = data.frame(
Model = c("SL-12", rep("", 9), "SL-12 (log)", rep("", 9),
"SL-12 (all)",rep("", 9), "SL-90", rep("", 9),
"SL-90 (log)", rep("", 9), "SL-90 (all)", rep("", 9),
"SAFE", rep("", 9)),
Decile = 1:10,
matrix(paste0(
round(OE_decile_CI_F2_df[,seq(1, ncol(OE_decile_CI_F2_df), by = 3)], 3), " (",
round(OE_decile_CI_F2_df[,seq(2, ncol(OE_decile_CI_F2_df), by = 3)], 3), ", ",
round(OE_decile_CI_F2_df[,seq(3, ncol(OE_decile_CI_F2_df), by = 3)], 3), ")"),
nrow = nrow(OE_decile_CI_F2_df), ncol = ncol(OE_decile_CI_F2_df)/3),
matrix(paste0(
round(OE_decile_CI_F3_df[,seq(1, ncol(OE_decile_CI_F3_df), by = 3)], 3), " (",
round(OE_decile_CI_F3_df[,seq(2, ncol(OE_decile_CI_F3_df), by = 3)], 3), ", ",
round(OE_decile_CI_F3_df[,seq(3, ncol(OE_decile_CI_F3_df), by = 3)], 3), ")"),
nrow = nrow(OE_decile_CI_F3_df), ncol = ncol(OE_decile_CI_F3_df)/3)
)
colnames(OE_decile_CI_table) =
c("", "Decile",
paste(c("FLINT", "NHANES-NAFLD"),
rep(c("(F2 or higher)", "(F3 or higher)"), each = 2)))
write.csv(OE_decile_CI_table, file = "figs_and_tabs/OE_decile_CI_table.csv")
knitr::kable(OE_decile_CI_table)| Decile | FLINT (F2 or higher) | NHANES-NAFLD (F2 or higher) | FLINT (F3 or higher) | NHANES-NAFLD (F3 or higher) | |
|---|---|---|---|---|---|
| SL-12 | 1 | 0.558 (0.134, 1.082) | 0.259 (0.028, 0.537) | 1.054 (0, 2.25) | 0.324 (0, 1.038) |
| 2 | 0.504 (0.11, 0.95) | 0.108 (0.018, 0.236) | 0.287 (0, 0.826) | 0.055 (0, 0.179) | |
| 3 | 1.03 (0.607, 1.427) | 0.418 (0.136, 0.772) | 1.007 (0.25, 1.787) | 0 (0, 0) | |
| 4 | 1.03 (0.686, 1.413) | 0.275 (0.098, 0.495) | 1.219 (0.585, 1.955) | 0.332 (0.05, 0.725) | |
| 5 | 1.074 (0.738, 1.404) | 0.34 (0.142, 0.608) | 1.123 (0.484, 1.809) | 0.229 (0.026, 0.545) | |
| 6 | 1.108 (0.762, 1.461) | 0.459 (0.253, 0.753) | 1.09 (0.481, 1.805) | 0.51 (0.092, 1.126) | |
| 7 | 1.192 (0.931, 1.417) | 0.441 (0.226, 0.72) | 0.792 (0.334, 1.257) | 0.129 (0.015, 0.301) | |
| 8 | 1.212 (1.021, 1.387) | 0.248 (0.123, 0.435) | 1.268 (0.839, 1.718) | 0.065 (0, 0.157) | |
| 9 | 1.015 (0.804, 1.195) | 0.623 (0.374, 0.859) | 0.892 (0.54, 1.243) | 0.925 (0.36, 1.49) | |
| 10 | 1.109 (0.995, 1.191) | 0.621 (0.456, 0.783) | 1.299 (1.089, 1.479) | 0.582 (0.35, 0.852) | |
| SL-12 (log) | 1 | 0.574 (0.142, 1.086) | 0.286 (0.031, 0.605) | 1.136 (0, 2.372) | 0.398 (0, 1.264) |
| 2 | 0.515 (0.109, 0.96) | 0.118 (0.02, 0.26) | 0.326 (0, 0.977) | 0.062 (0, 0.204) | |
| 3 | 1.019 (0.602, 1.423) | 0.445 (0.145, 0.832) | 1.123 (0.287, 2.003) | 0 (0, 0) | |
| 4 | 1.031 (0.682, 1.404) | 0.309 (0.11, 0.551) | 1.303 (0.62, 2.117) | 0.379 (0.059, 0.83) | |
| 5 | 1.055 (0.727, 1.387) | 0.391 (0.165, 0.688) | 1.244 (0.534, 2.007) | 0.256 (0.029, 0.608) | |
| 6 | 1.107 (0.765, 1.454) | 0.511 (0.288, 0.818) | 1.258 (0.561, 2.058) | 0.569 (0.101, 1.236) | |
| 7 | 1.18 (0.91, 1.396) | 0.499 (0.257, 0.809) | 0.881 (0.381, 1.379) | 0.148 (0.017, 0.34) | |
| 8 | 1.218 (1.022, 1.398) | 0.282 (0.14, 0.487) | 1.458 (0.954, 1.977) | 0.072 (0, 0.173) | |
| 9 | 1.031 (0.812, 1.216) | 0.715 (0.432, 0.985) | 0.994 (0.602, 1.405) | 1.05 (0.405, 1.682) | |
| 10 | 1.109 (0.992, 1.193) | 0.674 (0.496, 0.844) | 1.453 (1.219, 1.661) | 0.654 (0.4, 0.95) | |
| SL-12 (all) | 1 | 0.544 (0.135, 1.013) | 0.24 (0.026, 0.506) | 1.092 (0, 2.29) | 0.347 (0, 1.102) |
| 2 | 0.504 (0.108, 0.939) | 0.104 (0.018, 0.229) | 0.299 (0, 0.866) | 0.06 (0, 0.2) | |
| 3 | 0.999 (0.592, 1.384) | 0.409 (0.134, 0.752) | 0.958 (0.236, 1.714) | 0 (0, 0) | |
| 4 | 1.051 (0.704, 1.434) | 0.281 (0.101, 0.504) | 1.203 (0.577, 1.914) | 0.376 (0.057, 0.821) | |
| 5 | 1.066 (0.732, 1.399) | 0.356 (0.15, 0.625) | 1.07 (0.469, 1.721) | 0.266 (0.03, 0.623) | |
| 6 | 1.128 (0.782, 1.478) | 0.482 (0.267, 0.782) | 1.097 (0.486, 1.794) | 0.565 (0.1, 1.235) | |
| 7 | 1.195 (0.929, 1.419) | 0.46 (0.232, 0.748) | 0.76 (0.325, 1.201) | 0.147 (0.017, 0.339) | |
| 8 | 1.232 (1.032, 1.413) | 0.26 (0.128, 0.457) | 1.219 (0.801, 1.652) | 0.071 (0, 0.171) | |
| 9 | 1.035 (0.812, 1.225) | 0.675 (0.406, 0.933) | 0.857 (0.519, 1.203) | 1.064 (0.412, 1.722) | |
| 10 | 1.117 (1.005, 1.197) | 0.646 (0.473, 0.813) | 1.255 (1.056, 1.419) | 0.599 (0.366, 0.866) | |
| SL-90 | 1 | 0.545 (0.137, 1.032) | 0.262 (0.029, 0.553) | 1.289 (0, 2.712) | 0.407 (0, 1.292) |
| 2 | 0.517 (0.11, 0.969) | 0.108 (0.018, 0.236) | 0.356 (0, 1.053) | 0.065 (0, 0.214) | |
| 3 | 1.034 (0.614, 1.438) | 0.409 (0.134, 0.753) | 1.201 (0.292, 2.12) | 0 (0, 0) | |
| 4 | 1.068 (0.705, 1.463) | 0.275 (0.098, 0.491) | 1.486 (0.703, 2.409) | 0.405 (0.062, 0.893) | |
| 5 | 1.103 (0.759, 1.429) | 0.352 (0.148, 0.615) | 1.349 (0.589, 2.163) | 0.277 (0.031, 0.654) | |
| 6 | 1.146 (0.79, 1.515) | 0.473 (0.264, 0.77) | 1.337 (0.583, 2.232) | 0.596 (0.107, 1.301) | |
| 7 | 1.2 (0.943, 1.423) | 0.454 (0.23, 0.737) | 0.905 (0.384, 1.411) | 0.151 (0.017, 0.352) | |
| 8 | 1.254 (1.052, 1.433) | 0.264 (0.129, 0.462) | 1.502 (1, 2.03) | 0.077 (0, 0.187) | |
| 9 | 1.035 (0.815, 1.218) | 0.672 (0.402, 0.928) | 1.038 (0.629, 1.461) | 1.112 (0.431, 1.795) | |
| 10 | 1.128 (1.013, 1.211) | 0.651 (0.476, 0.821) | 1.484 (1.249, 1.678) | 0.689 (0.418, 1.006) | |
| SL-90 (log) | 1 | 0.563 (0.139, 1.048) | 0.254 (0.028, 0.533) | 1.085 (0, 2.259) | 0.292 (0, 0.919) |
| 2 | 0.5 (0.106, 0.941) | 0.11 (0.018, 0.242) | 0.335 (0, 0.98) | 0.055 (0, 0.181) | |
| 3 | 1 (0.586, 1.399) | 0.419 (0.135, 0.775) | 1.278 (0.312, 2.262) | 0 (0, 0) | |
| 4 | 1.01 (0.665, 1.379) | 0.287 (0.102, 0.515) | 1.591 (0.763, 2.572) | 0.382 (0.058, 0.842) | |
| 5 | 1.054 (0.731, 1.376) | 0.362 (0.153, 0.644) | 1.422 (0.623, 2.248) | 0.27 (0.031, 0.624) | |
| 6 | 1.1 (0.759, 1.443) | 0.494 (0.275, 0.803) | 1.404 (0.627, 2.375) | 0.672 (0.116, 1.558) | |
| 7 | 1.156 (0.892, 1.376) | 0.478 (0.243, 0.785) | 1.03 (0.428, 1.613) | 0.173 (0.02, 0.409) | |
| 8 | 1.216 (1.024, 1.397) | 0.265 (0.13, 0.461) | 1.62 (1.072, 2.171) | 0.091 (0, 0.218) | |
| 9 | 1.016 (0.803, 1.197) | 0.685 (0.413, 0.945) | 1.15 (0.697, 1.603) | 1.332 (0.516, 2.139) | |
| 10 | 1.104 (0.994, 1.184) | 0.657 (0.48, 0.827) | 1.659 (1.4, 1.881) | 0.799 (0.481, 1.17) | |
| SL-90 (all) | 1 | 0.545 (0.135, 1.037) | 0.224 (0.025, 0.47) | 1.018 (0, 2.091) | 0.312 (0, 0.991) |
| 2 | 0.504 (0.108, 0.946) | 0.099 (0.017, 0.218) | 0.317 (0, 0.919) | 0.056 (0, 0.186) | |
| 3 | 1.028 (0.609, 1.444) | 0.388 (0.127, 0.71) | 1.076 (0.259, 1.904) | 0 (0, 0) | |
| 4 | 1.061 (0.701, 1.443) | 0.268 (0.094, 0.484) | 1.308 (0.623, 2.111) | 0.365 (0.056, 0.797) | |
| 5 | 1.099 (0.762, 1.42) | 0.337 (0.141, 0.604) | 1.153 (0.502, 1.845) | 0.262 (0.029, 0.613) | |
| 6 | 1.16 (0.803, 1.52) | 0.465 (0.254, 0.772) | 1.186 (0.527, 2.001) | 0.576 (0.1, 1.3) | |
| 7 | 1.197 (0.936, 1.427) | 0.434 (0.22, 0.711) | 0.837 (0.354, 1.312) | 0.144 (0.016, 0.342) | |
| 8 | 1.247 (1.049, 1.433) | 0.247 (0.121, 0.43) | 1.309 (0.866, 1.761) | 0.075 (0, 0.178) | |
| 9 | 1.047 (0.824, 1.234) | 0.646 (0.387, 0.894) | 0.896 (0.544, 1.251) | 1.087 (0.412, 1.754) | |
| 10 | 1.134 (1.03, 1.214) | 0.633 (0.461, 0.798) | 1.297 (1.088, 1.473) | 0.619 (0.371, 0.904) | |
| SAFE | 1 | 0.957 (0.227, 1.886) | 0.596 (0.066, 1.231) | 0.718 (0, 1.507) | 0.235 (0, 0.732) |
| 2 | 0.652 (0.133, 1.201) | 0.171 (0.029, 0.378) | 0.13 (0, 0.398) | 0.028 (0, 0.095) | |
| 3 | 1.126 (0.658, 1.59) | 0.545 (0.179, 1.032) | 0.469 (0.095, 0.86) | 0 (0, 0) | |
| 4 | 1.164 (0.782, 1.555) | 0.344 (0.122, 0.616) | 0.621 (0.309, 0.998) | 0.145 (0.022, 0.327) | |
| 5 | 1.065 (0.727, 1.398) | 0.411 (0.172, 0.731) | 0.533 (0.201, 0.861) | 0.101 (0.011, 0.23) | |
| 6 | 1.014 (0.712, 1.315) | 0.499 (0.273, 0.823) | 0.477 (0.237, 0.731) | 0.227 (0.041, 0.536) | |
| 7 | 1.141 (0.874, 1.354) | 0.473 (0.25, 0.758) | 0.38 (0.162, 0.603) | 0.06 (0.007, 0.138) | |
| 8 | 1.144 (0.949, 1.302) | 0.255 (0.125, 0.435) | 0.746 (0.497, 0.998) | 0.029 (0, 0.068) | |
| 9 | 0.948 (0.731, 1.127) | 0.588 (0.361, 0.804) | 0.542 (0.317, 0.764) | 0.396 (0.157, 0.617) | |
| 10 | 1.026 (0.91, 1.112) | 0.562 (0.404, 0.71) | 0.944 (0.786, 1.071) | 0.31 (0.182, 0.454) |
We plot calibration by SAFE decile for the superlearner with 12 base models and the SAFE probability applied to the F2 or higher significant fibrosis outcome.
# FLINT F2
p1 = data.frame(
Model = factor(rep(c("Superlearner-F2", "SAFE"), each = 10),
levels = c("Superlearner-F2", "SAFE")),
Decile = 1:10,
x = rep(1:10, 2) + rep(c(-0.15, 0.15), each = 10),
data.frame(
OE_decile_CI_F2_df[rep(1:10, 2) + rep(c(20, 60), each = 10), 1:3]) %>%
setNames(c("est", "lo", "hi"))
) %>%
ggplot(aes(x = x, y = est, color = Model)) +
geom_pointrange(aes(ymin = lo, ymax = hi), size = 0.5, stroke = 0.5) +
geom_hline(yintercept = 1, color = "grey") +
scale_x_continuous(breaks = 1:10) +
xlab("Decile") + ylab("Calibration (Observed / Expected)") +
ggtitle("FLINT calibration (F2 or higher)")
ggsave("figs_and_tabs/OE_decile_flint_F2.png", p1, width = 6.5, height = 4)
p1
# NHANES F2
p1 = data.frame(
Model = factor(rep(c("Superlearner-F2", "SAFE"), each = 10),
levels = c("Superlearner-F2", "SAFE")),
Decile = 1:10,
x = rep(1:10, 2) + rep(c(-0.15, 0.15), each = 10),
data.frame(
OE_decile_CI_F2_df[rep(1:10, 2) + rep(c(20, 60), each = 10), 4:6]) %>%
setNames(c("est", "lo", "hi"))
) %>%
ggplot(aes(x = x, y = est, color = Model)) +
geom_pointrange(aes(ymin = lo, ymax = hi), size = 0.5, stroke = 0.5) +
geom_hline(yintercept = 1, color = "grey") +
scale_x_continuous(breaks = 1:10) +
xlab("Decile") + ylab("Calibration (Observed / Expected)") +
ggtitle("NHANES-NAFLD calibration (F2 or higher)")
ggsave("figs_and_tabs/OE_decile_nhanes_F2.png", p1, width = 6.5, height = 4)
p1
Spiegelhalter z test as a proxy for calibration. Significance (i.e. improper calibration) is determined at an alpha level of 0.05. Note that the calculation for NHANES-NAFLD cohort is not weighted, but the complex survey design should probably be accounted for.
# Function to calculate Spiegelhalter z test for multiple scores
calc_Spiegelhalter_z = function(preds, response, alpha = 0.05) {
sapply(preds, function(x) {
z = sum((response - x) * (1-2*x)) / sqrt(sum((1-2*x)^2 * x * (1-x)))
sig = abs(z) > qnorm(1 - alpha/2)
return(data.frame(z = z, sig = sig))
})
}
# FLINT
calc_Spiegelhalter_z(scores_flint[,prob_F2_idx], flintY2)## SL_12_F2 SL_12_log_F2 SL_12_all_F2 SL_90_F2 SL_90_log_F2 SL_90_all_F2
## z -2.298313 -2.548487 -2.796165 -2.950836 -2.035195 -2.999257
## sig TRUE TRUE TRUE TRUE TRUE TRUE
## SAFE_prob
## z -1.126894
## sig FALSEcalc_Spiegelhalter_z(scores_flint[,prob_F3_idx], flintY3)## SL_12_F3 SL_12_log_F3 SL_12_all_F3 SL_90_F3 SL_90_log_F3 SL_90_all_F3
## z -0.4616254 0.4006335 -0.6981216 0.9122109 1.519016 -0.07791858
## sig FALSE FALSE FALSE FALSE FALSE FALSE
## SAFE_prob
## z 4.138205
## sig TRUE# NHANES
calc_Spiegelhalter_z(scores_nhanes[,prob_F2_idx], nhanesY2)## SL_12_F2 SL_12_log_F2 SL_12_all_F2 SL_90_F2 SL_90_log_F2 SL_90_all_F2
## z -4.19961 -9.008777 -6.944887 -7.454138 -6.922753 -4.503174
## sig TRUE TRUE TRUE TRUE TRUE TRUE
## SAFE_prob
## z -2.53807
## sig TRUEcalc_Spiegelhalter_z(scores_nhanes[,prob_F3_idx], nhanesY3)## SL_12_F3 SL_12_log_F3 SL_12_all_F3 SL_90_F3 SL_90_log_F3 SL_90_all_F3
## z -9.798676 -9.127421 -8.801286 -8.821021 -8.558064 -8.990836
## sig TRUE TRUE TRUE TRUE TRUE TRUE
## SAFE_prob
## z -3.888703
## sig TRUECox intercept and slope method, with interpretation as follows:
- Slope = direction of miscalibration, where 1 is perfect, greater than 1 means underestimating high risk/overestimating low risk, and lower than 1 means overestimating high risk/underestimating low risk
- Intercept = overall miscalibration, where 0 is good calibration, greater than 0 means average understimation, and lower than 0 means average overestimation.
We used generalized linear models for complex survey designs for NHANES-NAFLD, but something still seems slightly odd about the results, especially the intercepts.
# Function to calculate Cox intercept and slope for multiple scores
cox_intercept_slope = function(preds, response) {
sapply(preds, function(x) {
my_df = data.frame(o = response, e = x)
my_df$logite = log(my_df$e / (1 - my_df$e))
fit = glm(o ~ I(logite), data = my_df, family = binomial(link = "logit"))
est = coef(fit)
ci = confint(fit)
return(c(Intercept = paste0(round(est[1], 3), " (", round(ci[1,1], 3),
", ", round(ci[1,2], 3), ")"),
Slope = paste0(round(est[2], 3), " (", round(ci[2,1], 3),
", ", round(ci[2,2], 3), ")")))
})
}
# Function to calculate Cox intercept and slope for multiple scores,
# while accounting for sampling weights
svy_cox_intercept_slope = function(preds, response, design) {
sapply(preds, function(x) {
my_df = data.frame(o = response, e = x)
my_df$logite = log(my_df$e / (1 - my_df$e))
design = svydesign(id = ~nhanes_clusters,
strata = ~nhanes_strata,
weights = ~nhanes_weights, nest=TRUE,
data = my_df)
fit = svyglm(o ~ I(logite), design = design, family = binomial(link = "logit"))
est = coef(fit)
ci = confint(fit)
return(c(Intercept = paste0(round(est[1], 3), " (", round(ci[1,1], 3),
", ", round(ci[1,2], 3), ")"),
Slope = paste0(round(est[2], 3), " (", round(ci[2,1], 3),
", ", round(ci[2,2], 3), ")")))
})
}
# FLINT
cox_intercept_slope(scores_flint[,prob_F2_idx], flintY2) %>% t() %>% knitr::kable()| Intercept | Slope | |
|---|---|---|
| SL_12_F2 | 0.006 (-0.287, 0.294) | 1.459 (1.077, 1.883) |
| SL_12_log_F2 | -0.006 (-0.303, 0.284) | 1.531 (1.134, 1.973) |
| SL_12_all_F2 | 0.002 (-0.294, 0.292) | 1.594 (1.186, 2.048) |
| SL_90_F2 | 0.06 (-0.231, 0.348) | 1.626 (1.21, 2.09) |
| SL_90_log_F2 | -0.023 (-0.32, 0.267) | 1.405 (1.035, 1.816) |
| SL_90_all_F2 | 0.058 (-0.235, 0.346) | 1.635 (1.222, 2.096) |
| SAFE_prob | 0.104 (-0.184, 0.389) | 1.158 (0.856, 1.493) |
cox_intercept_slope(scores_flint[,prob_F3_idx], flintY3) %>% t() %>% knitr::kable()| Intercept | Slope | |
|---|---|---|
| SL_12_F3 | 0.253 (-0.109, 0.63) | 1.194 (0.868, 1.556) |
| SL_12_log_F3 | 0.567 (0.151, 1.005) | 1.3 (0.95, 1.689) |
| SL_12_all_F3 | 0.191 (-0.16, 0.557) | 1.188 (0.87, 1.54) |
| SL_90_F3 | 0.689 (0.252, 1.152) | 1.314 (0.963, 1.703) |
| SL_90_log_F3 | 0.976 (0.481, 1.503) | 1.451 (1.07, 1.871) |
| SL_90_all_F3 | 0.31 (-0.062, 0.701) | 1.173 (0.859, 1.521) |
| SAFE_prob | -1.17 (-1.523, -0.845) | 0.956 (0.67, 1.274) |
# NHANES
svy_cox_intercept_slope(scores_nhanes[,prob_F2_idx], nhanesY2, nhanes_design) %>% t() %>% knitr::kable()| Intercept | Slope | |
|---|---|---|
| SL_12_F2 | -1.348 (-1.587, -1.109) | 0.946 (0.579, 1.312) |
| SL_12_log_F2 | -1.014 (-1.268, -0.76) | 1.246 (0.851, 1.641) |
| SL_12_all_F2 | -1.245 (-1.497, -0.992) | 1.047 (0.676, 1.418) |
| SL_90_F2 | -1.225 (-1.462, -0.987) | 1.087 (0.715, 1.459) |
| SL_90_log_F2 | -1.207 (-1.448, -0.966) | 1.028 (0.635, 1.421) |
| SL_90_all_F2 | -1.377 (-1.61, -1.143) | 0.874 (0.485, 1.264) |
| SAFE_prob | -1.262 (-1.537, -0.986) | 0.85 (0.591, 1.109) |
svy_cox_intercept_slope(scores_nhanes[,prob_F3_idx], nhanesY3, nhanes_design) %>% t() %>% knitr::kable()| Intercept | Slope | |
|---|---|---|
| SL_12_F3 | -1.048 (-1.61, -0.485) | 1.12 (0.797, 1.442) |
| SL_12_log_F3 | -0.839 (-1.461, -0.217) | 1.138 (0.792, 1.483) |
| SL_12_all_F3 | -0.917 (-1.545, -0.29) | 1.122 (0.797, 1.447) |
| SL_90_F3 | -0.618 (-1.277, 0.04) | 1.23 (0.874, 1.586) |
| SL_90_log_F3 | -0.441 (-1.181, 0.3) | 1.254 (0.85, 1.658) |
| SL_90_all_F3 | -0.94 (-1.547, -0.333) | 1.094 (0.777, 1.411) |
| SAFE_prob | -2.422 (-2.844, -2.001) | 1.095 (0.664, 1.526) |
SAFE sensitivity, specificity, PPV, and NPV
We calculate the sensitivity, specificity, PPV, and NPV for SAFE at the thresholds of 100 (SAFE >= 100 corresponds to high risk) and 0 (SAFE < 0 corresponds to low risk).
# FLINT
data.frame(
Outcome = rep(c("F2 or higher", "F3 or higher"), each = 2),
rbind(coords(roc(flintY2, scores_flint[,"SAFE"]),
c(0, 100), input = "threshold",
ret = c("threshold", "sensitivity","specificity","ppv","npv")),
coords(roc(flintY3, scores_flint[,"SAFE"]),
c(0, 100), input = "threshold",
ret = c("threshold", "sensitivity","specificity","ppv","npv")))
) %>%
setNames(c("Outcome", "SAFE Threshold",
"Sensitivity", "Specificity", "PPV", "NPV")) %>%
mutate_if(is.numeric, round, 3) %>% knitr::kable()| Outcome | SAFE Threshold | Sensitivity | Specificity | PPV | NPV |
|---|---|---|---|---|---|
| F2 or higher | 0 | 0.962 | 0.261 | 0.651 | 0.829 |
| F2 or higher | 100 | 0.780 | 0.649 | 0.761 | 0.673 |
| F3 or higher | 0 | 0.967 | 0.178 | 0.370 | 0.914 |
| F3 or higher | 100 | 0.822 | 0.506 | 0.454 | 0.850 |
# NHANES-NAFLD
data.frame(
Outcome = rep(c("F2 or higher", "F3 or higher"), each = 2),
rbind(coords(roc(nhanesY2, scores_nhanes[,"SAFE"]),
c(0, 100), input = "threshold",
ret = c("threshold", "sensitivity","specificity","ppv","npv")),
coords(roc(nhanesY3, scores_nhanes[,"SAFE"]),
c(0, 100), input = "threshold",
ret = c("threshold", "sensitivity","specificity","ppv","npv")))
) %>%
setNames(c("Outcome", "SAFE Threshold",
"Sensitivity", "Specificity", "PPV", "NPV")) %>%
mutate_if(is.numeric, round, 3) %>% knitr::kable()| Outcome | SAFE Threshold | Sensitivity | Specificity | PPV | NPV |
|---|---|---|---|---|---|
| F2 or higher | 0 | 0.842 | 0.455 | 0.224 | 0.939 |
| F2 or higher | 100 | 0.480 | 0.822 | 0.335 | 0.894 |
| F3 or higher | 0 | 0.910 | 0.427 | 0.083 | 0.988 |
| F3 or higher | 100 | 0.657 | 0.799 | 0.157 | 0.976 |
Scatterplots
This is a plot of the superlearner (90 base models, log-transformed predictors) vs. SAFE scores for each cohort. Those with significant fibrosis are colored pink and those without are blue. The Spearman’s correlation coefficient is reported in the lower right; the NHANES-NAFLD correlation accounts for the survey design.
# Function to calculate Spearman's correlation for complex surveys
svyspear = function(x, y, design) {
vmat = as.matrix(svyvar(~rank(x)+rank(y), design))
attr(vmat, "var") = NULL #remove extra info
attr(vmat, "statistic") = NULL #remove extra info
cov2cor(vmat)[1,2]
}
# Combine scores and outcomes for all cohorts
stacked_scores_df = rbind(
data.frame(setNames(scores_flint[,c("SL_12_all_F2", "SAFE")], c("SL", "SAFE")),
sigfibro = flint$sigfibro2,
outcome = "Superlearner (F2 or higher)",
cohort = "FLINT"),
data.frame(setNames(scores_flint[,c("SL_12_all_F3", "SAFE")], c("SL", "SAFE")),
sigfibro = flint$sigfibro3,
outcome = "Superlearner (F3 or higher)",
cohort = "FLINT"),
data.frame(setNames(scores_nhanes[,c("SL_12_all_F2", "SAFE")], c("SL", "SAFE")),
sigfibro = nhanes$sigfibro2,
outcome = "Superlearner (F2 or higher)",
cohort = "NHANES-NAFLD"),
data.frame(setNames(scores_nhanes[,c("SL_12_all_F3", "SAFE")], c("SL", "SAFE")),
sigfibro = nhanes$sigfibro3,
outcome = "Superlearner (F3 or higher)",
cohort = "NHANES-NAFLD")
)
# Spearman's correlation between superlearner and SAFE
sl_safe_spear = rbind(
data.frame(
outcome = "Superlearner (F2 or higher)",
cohort = c("FLINT", "NHANES-NAFLD"),
Cor = c(with(scores_flint, cor(SL_12_all_F2, SAFE)),
with(scores_nhanes, svyspear(SL_12_all_F2, SAFE, nhanes_design)))
),
data.frame(
outcome = "Superlearner (F3 or higher)",
cohort = c("FLINT", "NHANES-NAFLD"),
Cor = c(with(scores_flint, cor(SL_12_all_F3, SAFE)),
with(scores_nhanes, svyspear(SL_12_all_F3, SAFE, nhanes_design)))
)
)
# Plot superlearner against SAFE scores
p1 = stacked_scores_df %>%
ggplot(aes(x = SAFE, y = SL, col = as.factor(sigfibro))) +
geom_point(alpha = 0.5) +
facet_grid(cohort~outcome) +
geom_label(data = sl_safe_spear, aes(label = round(Cor, 2), col = "black"),
x = Inf, y = -Inf, vjust = -0.2, hjust = "inward",
alpha = 0.8, size = 3, show.legend = FALSE) +
scale_color_discrete(name = "Significant Fibrosis",
limits = c("1", "0"), labels = c("Yes", "No")) +
ylab("Superlearner") + theme(legend.position="bottom")
ggsave("figs_and_tabs/sl_vs_safe.png", width = 6, height = 6.5)
p1
We plot log-transformed LSM against the superlearner (90 base models, log-transformed predictors), APRI, BARD, FIB4, Forns, NFS, and SAFE scores in the NHANES-NAFLD cohort.
# LSM for NHANES-NAFLD
nhanes_lsm_df = data.frame(
LSM = nhanes$LSM,
scores_nhanes[,c(mod_sub_names_F2[1], mod_sub_names_F3)],
Weights = nhanes_weights) %>%
rename("FIB-4" = "FIB4",
"Superlearner\n(F2 or higher)" = "SL_12_all_F2",
"Superlearner\n(F3 or higher)" = "SL_12_all_F3") %>%
melt(id.vars = c("LSM", "Weights"), variable.name = "Model",
value.name = "Score")
# Spearman's correlation with LSM
nhanes_lsm_spear = nhanes_lsm_df %>% group_by(Model) %>%
summarise("NHANES-NAFLD" = svyspear(LSM, Score, nhanes_design))
# Plot LSM against scores
p1 = nhanes_lsm_df %>%
ggplot(aes(x = Score, y = LSM, col = Model)) +
geom_point(alpha = 0.5) +
geom_smooth(method = 'loess', se = FALSE, lwd = 0.5) +
scale_y_log10() +
geom_label(data = nhanes_lsm_spear, aes(label = round(`NHANES-NAFLD`, 2)),
x = Inf, y = -Inf, vjust = -0.2, hjust = "inward",
alpha = 0.8, size = 3) +
facet_wrap(~Model, scales = "free_x", ncol = 4) +
ylab("Liver Stiffness Measure (kPa)") +
guides(col = "none")
ggsave("figs_and_tabs/nhanes_lsm.png", p1, width = 7, height = 4.5)
p1
We calculate the FAST, Agile3+ and Agile4 scores for the NHANES-NAFLD cohort.
# Expit function
expit = function(x) {
exp(x) / (1 + exp(x))
}
# NHANES-NAFLD FAST, Agile3+, and Agile4 scores
nhanes = nhanes %>%
mutate(FAST = expit(-1.65 + 1.07*log(LSM) +
2.66*10^(-8)*CAP^3 - 63.3/AST),
Agile3 = -3.92368 + 2.29714 * log(LSM) - 0.00902 * PLAT -
0.98633 / (AST/ALT) + 1.08636 * DIAB2 - 0.38581 * MALE + 0.03018 * AGE,
Agile4 = 7.50139 - 15.42498 / sqrt(LSM) - 0.01378 * PLAT -
1.41149 / (AST/ALT) - 0.53281 * MALE + 0.41741 * DIAB2)This figure plots the log-transformed Fibro-AST score against the superlearner (90 base models, log-transformed predictors), APRI, BARD, FIB4, Forns, NFS, and SAFE scores.
# Fibro-AST for NHANES-NAFLD
nhanes_fast_df = data.frame(FAST = nhanes$FAST,
scores_nhanes[,c(mod_sub_names_F2[1], mod_sub_names_F3)],
Weights = nhanes_weights) %>%
rename("FIB-4" = "FIB4",
"Superlearner\n(F2 or higher)" = "SL_12_all_F2",
"Superlearner\n(F3 or higher)" = "SL_12_all_F3") %>%
melt(id.vars = c("FAST", "Weights"), variable.name = "Model",
value.name = "Score")
# Spearman's correlation with FAST
nhanes_fast_spear = nhanes_fast_df %>% group_by(Model) %>%
summarise("NHANES-NAFLD" = svyspear(FAST, Score, nhanes_design))
# Plot FAST against other scores
p1 = nhanes_fast_df %>%
ggplot(aes(x = Score, y = FAST, col = Model)) +
geom_point(alpha = 0.5) +
geom_smooth(method = 'loess', se = FALSE, lwd = 0.5) +
scale_y_log10() +
geom_label(data = nhanes_fast_spear, aes(label = round(`NHANES-NAFLD`, 2)),
x = Inf, y = -Inf, vjust = -0.2, hjust = "inward",
alpha = 0.8, size = 3) +
facet_wrap(~Model, scales = "free_x", ncol = 4) +
ylab("FibroScan-AST") +
guides(col = "none")
ggsave("figs_and_tabs/nhanes_fast.png", p1, width = 7, height = 4.5)
p1
We plot Agile3+ against the superlearner (90 base models, log-transformed predictors), BARD, APRI, FIB4, Forns, NFS, and SAFE scores.
# Agile3+ for NHANES-NAFLD
nhanes_agile3_df = data.frame(Agile3 = nhanes$Agile3,
scores_nhanes[,c(mod_sub_names_F2[1], mod_sub_names_F3)],
Weights = nhanes_weights) %>%
rename("FIB-4" = "FIB4",
"Superlearner\n(F2 or higher)" = "SL_12_all_F2",
"Superlearner\n(F3 or higher)" = "SL_12_all_F3") %>%
melt(id.vars = c("Agile3", "Weights"), variable.name = "Model",
value.name = "Score")
# Spearman's correlation with Agile3+
nhanes_agile3_spear = nhanes_agile3_df %>% group_by(Model) %>%
summarise("NHANES-NAFLD" = svyspear(Agile3, Score, nhanes_design))
# Plot Agile3+ against other scores
p1 = nhanes_agile3_df %>%
ggplot(aes(x = Score, y = Agile3, col = Model)) +
geom_point(alpha = 0.5) +
geom_smooth(method = 'loess', se = FALSE, lwd = 0.5) +
geom_label(data = nhanes_agile3_spear, aes(label = round(`NHANES-NAFLD`, 2)),
x = Inf, y = -Inf, vjust = -0.2, hjust = "inward",
alpha = 0.8, size = 3) +
facet_wrap(~Model, scales = "free_x", ncol = 4) +
ylab("Agile3+") +
guides(col = "none")
ggsave("figs_and_tabs/nhanes_agile3.png", p1, width = 7, height = 4.5)
p1
We plot Agile4 against the superlearner (90 base models, log-transformed predictors), APRI, BARD, FIB4, Forns, NFS, and SAFE scores.
# Agile4 for NHANES-NAFLD
nhanes_agile4_df = data.frame(Agile4 = nhanes$Agile4,
scores_nhanes[,c(mod_sub_names_F2[1], mod_sub_names_F3)],
Weights = nhanes_weights) %>%
rename("FIB-4" = "FIB4",
"Superlearner\n(F2 or higher)" = "SL_12_all_F2",
"Superlearner\n(F3 or higher)" = "SL_12_all_F3") %>%
melt(id.vars = c("Agile4", "Weights"), variable.name = "Model",
value.name = "Score")
# Spearman's correlation with Agile4
nhanes_agile4_spear = nhanes_agile4_df %>% group_by(Model) %>%
summarise("NHANES-NAFLD" = svyspear(Agile4, Score, nhanes_design))
# Plot Agile4 against other scores
p1 = nhanes_agile4_df %>%
ggplot(aes(x = Score, y = Agile4, col = Model)) +
geom_point(alpha = 0.5) +
geom_smooth(method = 'loess', se = FALSE, lwd = 0.5) +
geom_label(data = nhanes_agile4_spear, aes(label = round(`NHANES-NAFLD`, 2)),
x = Inf, y = -Inf, vjust = -0.2, hjust = "inward",
alpha = 0.8, size = 3) +
facet_wrap(~Model, scales = "free_x", ncol = 4) +
ylab("Agile4") +
guides(col = "none")
ggsave("figs_and_tabs/nhanes_agile4.png", p1, width = 7, height = 4.5)
p1
Standardized means by quartile
First, for each validation cohort, we standardize the continuous variables by subtracting the mean and dividing by the standard deviation. In the NHANES-NAFLD cohort, we use means and standard deviations calculated for complex survey designs, which incorporate sampling weights, clusters, and strata. Then, for each quartile defined by the superlearner and SAFE scores (weighted quartiles for NHANES-NAFLD), we can calculate the standardized mean of each continuous predictor and the mean of each binary predictor. Confidence intervals were obtained using the usual normal approximation for continuous covariates; we used Wald confidence intervals for the binary covariates. Again, the NHANES-NAFLD mean and confidence interval reported are based on the appropriate estimates for complex survey designs.
# Functions to calculate standard errors for continuous and binary variables
SE = function(x) {
sd(x) / sqrt(length(x))
}
SE_binomial = function(x) {
sqrt(mean(x) * (1 - mean(x)) / length(x))
}
# Function to calculate weighted standard errors
svySE = function(x, design) {
data.frame(svymean(x, design))$SE
}
# FLINT
# Standardize all continuous variables
std_flint_df = sapply(1:length(pred_vars), function(i) {
if (!(pred_vars[i] %in% pred_vars_binary)) {
return(scale(flint[,pred_vars[i]]))
} else {
return(flint[,pred_vars[i]])
}
})
colnames(std_flint_df) = pred_vars
# Put together data frame with all predictors, SAFE, and superlearner quartiles
std_flint_df = data.frame(
std_flint_df,
quantiles_SL_F2 = cut(scores_flint$SL_12_all_F2,
breaks = c(-Inf, quantile(scores_flint$SL_12_all_F2,
c(0.25, 0.5, 0.75)), Inf),
labels = 1:4),
quantiles_SL_F3 = cut(scores_flint$SL_12_all_F3,
breaks = c(-Inf, quantile(scores_flint$SL_12_all_F3,
c(0.25, 0.5, 0.75)), Inf),
labels = 1:4),
quantiles_SAFE = cut(scores_flint$SAFE,
breaks = c(-Inf, quantile(scores_flint$SAFE,
c(0.25, 0.5, 0.75)), Inf),
labels = 1:4))
# Summarize data by quartile
std_flint_plot_df = rbind(
data.frame(
std_flint_df %>% rename(Quantile = "quantiles_SL_F2") %>%
group_by(Quantile) %>%
summarise(across(-c(quantiles_SAFE, quantiles_SL_F3),
list(est = ~mean(.x),
lo = ~(mean(.x) - qnorm(.975) *
ifelse(length(unique(.x)) == 2 &&
all(sort(unique(.x)) == c(0, 1)),
SE_binomial(.x), SE(.x))),
hi = ~(mean(.x) + qnorm(.975) *
ifelse(length(unique(.x)) == 2 &&
all(sort(unique(.x)) == c(0, 1)),
SE_binomial(.x), SE(.x)))))) %>%
pivot_longer(cols = -Quantile,
names_sep = "_",
names_to = c("variable", ".value")),
Score = "Superlearner\n(F2 or higher)"),
data.frame(
std_flint_df %>% rename(Quantile = "quantiles_SL_F3") %>%
group_by(Quantile) %>%
summarise(across(-c(quantiles_SAFE, quantiles_SL_F2),
list(est = ~mean(.x),
lo = ~(mean(.x) - qnorm(.975) *
ifelse(length(unique(.x)) == 2 &&
all(sort(unique(.x)) == c(0, 1)),
SE_binomial(.x), SE(.x))),
hi = ~(mean(.x) + qnorm(.975) *
ifelse(length(unique(.x)) == 2 &&
all(sort(unique(.x)) == c(0, 1)),
SE_binomial(.x), SE(.x)))))) %>%
pivot_longer(cols = -Quantile,
names_sep = "_",
names_to = c("variable", ".value")),
Score = "Superlearner\n(F3 or higher)"),
data.frame(
std_flint_df %>% rename(Quantile = "quantiles_SAFE") %>%
group_by(Quantile) %>%
summarise(across(-c(quantiles_SL_F2, quantiles_SL_F3),
list(est = ~mean(.x),
lo = ~(mean(.x) - qnorm(.975) *
ifelse(length(unique(.x)) == 2 &&
all(sort(unique(.x)) == c(0, 1)),
SE_binomial(.x), SE(.x))),
hi = ~(mean(.x) + qnorm(.975) *
ifelse(length(unique(.x)) == 2 &&
all(sort(unique(.x)) == c(0, 1)),
SE_binomial(.x), SE(.x)))))) %>%
pivot_longer(cols = -Quantile,
names_sep = "_",
names_to = c("variable", ".value")),
Score = "SAFE")
)
# NHANES-NAFLD
# Standardize all continuous variables
std_nhanes_df = sapply(1:length(pred_vars), function(i) {
if (!(pred_vars[i] %in% pred_vars_binary)) {
return((nhanes[,pred_vars[i]] - weighted.mean(nhanes[,pred_vars[i]], nhanes_weights)) /
sqrt(svyvar(nhanes[,pred_vars[i]], nhanes_design)))
} else {
return(nhanes[,pred_vars][,i])
}
})
colnames(std_nhanes_df) = pred_vars
# Put together data frame with all predictors, SAFE, and superlearner quartiles
std_nhanes_df = data.frame(
std_nhanes_df,
quantiles_SL_F2 = cut(scores_nhanes$SL_12_all_F2,
breaks = c(-Inf, Quantile(scores_nhanes$SL_12_all_F2,
nhanes_weights,
c(0.25, 0.5, 0.75)), Inf),
labels = 1:4),
quantiles_SL_F3 = cut(scores_nhanes$SL_12_all_F3,
breaks = c(-Inf, Quantile(scores_nhanes$SL_12_all_F3,
nhanes_weights,
c(0.25, 0.5, 0.75)), Inf),
labels = 1:4),
quantiles_SAFE = cut(scores_nhanes$SAFE,
breaks = c(-Inf, Quantile(scores_nhanes$SAFE,
nhanes_weights,
c(0.25, 0.5, 0.75)), Inf),
labels = 1:4),
weights = nhanes_weights)
# Subset complex survey designs by quartile
nhanes_design_by_quantiles_SL_F2 = lapply(1:4, function(i) {
subset(nhanes_design, std_nhanes_df$quantiles_SL_F2 == i)
})
nhanes_design_by_quantiles_SL_F3 = lapply(1:4, function(i) {
subset(nhanes_design, std_nhanes_df$quantiles_SL_F3 == i)
})
nhanes_design_by_quantiles_SAFE = lapply(1:4, function(i) {
subset(nhanes_design, std_nhanes_df$quantiles_SAFE == i)
})
# Summarize data by quartile
std_nhanes_plot_df = rbind(
data.frame(
std_nhanes_df %>% rename(Quantile = "quantiles_SL_F2") %>%
group_by(Quantile) %>%
summarise(across(
-c(weights, quantiles_SAFE, quantiles_SL_F3),
list(est = ~weighted.mean(.x, weights),
lo = ~(weighted.mean(.x, weights) - qnorm(.975) *
svySE(.x, nhanes_design_by_quantiles_SL_F2[[unique(Quantile)]])),
hi = ~(weighted.mean(.x, weights) + qnorm(.975) *
svySE(.x, nhanes_design_by_quantiles_SL_F2[[unique(Quantile)]]))))) %>%
pivot_longer(cols = -Quantile,
names_sep = "_",
names_to = c("variable", ".value")),
Score = "Superlearner\n(F2 or higher)"),
data.frame(
std_nhanes_df %>% rename(Quantile = "quantiles_SL_F3") %>%
group_by(Quantile) %>%
summarise(across(
-c(weights, quantiles_SAFE, quantiles_SL_F2),
list(est = ~weighted.mean(.x, weights),
lo = ~(weighted.mean(.x, weights) - qnorm(.975) *
svySE(.x, nhanes_design_by_quantiles_SL_F3[[unique(Quantile)]])),
hi = ~(weighted.mean(.x, weights) + qnorm(.975) *
svySE(.x, nhanes_design_by_quantiles_SL_F3[[unique(Quantile)]]))))) %>%
pivot_longer(cols = -Quantile,
names_sep = "_",
names_to = c("variable", ".value")),
Score = "Superlearner\n(F3 or higher)"),
data.frame(
std_nhanes_df %>% rename(Quantile = "quantiles_SAFE") %>%
group_by(Quantile) %>%
summarise(across(
-c(weights, quantiles_SL_F2, quantiles_SL_F3),
list(est = ~weighted.mean(.x, weights),
lo = ~(weighted.mean(.x, weights) - qnorm(.975) *
svySE(.x, nhanes_design_by_quantiles_SAFE[[unique(Quantile)]])),
hi = ~(weighted.mean(.x, weights) + qnorm(.975) *
svySE(.x, nhanes_design_by_quantiles_SAFE[[unique(Quantile)]]))))) %>%
pivot_longer(cols = -Quantile,
names_sep = "_",
names_to = c("variable", ".value")),
Score = "SAFE")
)
# Plot standardized means for each quartile, cohort, and score
p1 = rbind(data.frame(std_flint_plot_df, Cohort = "FLINT"),
data.frame(std_nhanes_plot_df, Cohort = "NHANES-NAFLD")) %>%
mutate(y = rep(rep(match(pred_vars, dat_dict$short), 4) +
rep(seq(-0.2, 0.2, length = 4), each = length(pred_vars)), 2*3)) %>%
ggplot(aes(x = as.numeric(est), y = y, color = Quantile)) +
geom_pointrange(aes(xmin = lo, xmax = hi), size = 0.25) +
facet_grid(Cohort ~
factor(Score, levels = c("Superlearner\n(F2 or higher)",
"Superlearner\n(F3 or higher)", "SAFE"))) +
scale_color_manual(values = c("navy", "dodgerblue", "mediumpurple","firebrick")) +
scale_y_continuous(breaks = match(pred_vars, dat_dict$short),
labels = dat_dict$long[match(pred_vars, dat_dict$short)],
expand = c(0.015, 0.015)) +
xlab("Standardized mean within quartile") + ylab("") +
theme(legend.position="bottom")
ggsave("figs_and_tabs/std_means.png", p1, width = 9, height = 8)
p1
# Plot standardized means for each quartile, cohort, and score for F2 only
p1 = rbind(data.frame(std_flint_plot_df, Cohort = "FLINT"),
data.frame(std_nhanes_plot_df, Cohort = "NHANES-NAFLD")) %>%
filter(Score != "Superlearner\n(F3 or higher)") %>%
mutate(y = rep(rep(match(pred_vars, dat_dict$short), 4) +
rep(seq(-0.2, 0.2, length = 4), each = length(pred_vars)), 2*2),
Score = fct_recode(Score, "Superlearner" = "Superlearner\n(F2 or higher)")) %>%
ggplot(aes(x = as.numeric(est), y = y, color = Quantile)) +
geom_pointrange(aes(xmin = lo, xmax = hi), size = 0.25) +
facet_grid(Cohort ~ Score) +
scale_color_grey() +
scale_y_continuous(breaks = match(pred_vars, dat_dict$short),
labels = dat_dict$long[match(pred_vars, dat_dict$short)],
expand = c(0.015, 0.015)) +
xlab("Standardized mean within quartile") + ylab("") +
theme(legend.position="bottom")
ggsave("figs_and_tabs/std_means_F2.png", p1, width = 6.5, height = 7)
p1